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1  Introduction 


*  .f 


Ultrasonic  diffraction  tomography  (UDT)  [1-3]  can  be  viewed  as  a  generalization  of  X-ray 
tomography  where  X-rays  have  been  replaced  with  an  acoustical  wavefield.  Because  UDT 
is  non-invasive,  free  of  radiation  hazard,  and  reproducible,  it  is  potentially  an  excellent  tool 
for  imaging  of  breast  cancer  [4,5].  While  UDT  promises  several  potentially  important  advan¬ 
tages  over  conventional  ultrasonic  imaging  and  has  found  important  uses  in  a  wide  variety  of 
engineering  and  scientific  disciplines,  its  application  to  imaging  of  breast  cancer  still  remains 
largely  unexplored.  The  broad  objective  of  the  proposed  project  is  to  investigate,  develop,  and 
evaluate  computationally  efficient  [6]  and  statistically  optimal  [7-9]  algorithms  for  accurate 
image  reconstruction  in  three-dimensional  (3D)  UDT  imaging  of  the  breast  cancer.  In  the  last 
year,  our  research  on  this  project  has,  we  believe,  been  successful  and  productive.  The  report 
below  summarizes  our  research  activities  and  results  on  the  project  to  date. 


2  Body 

Our  research  activities  on  the  project  to  date  can  be  grouped  naturally  into  4  categories.  The 
first  was  the  investigation  of  efficient  linear  algorithms  for  image  reconstruction  from  3D  data 
and  from  minimum  scan  data.  The  second  was  the  development  of  efficient  nonlinear  algo¬ 
rithms  for  2D  and  3D  image  reconstructions.  The  third  was  the  applications  of  the  developed 
algorithms  to  simulated  and  experimental  data  for  evaluation  of  their  performance.  Finally, 
as  a  by  product,  we  developed  short-scan  reconstruction  algorithms  for  reflection-mode  ultra¬ 
sound  tomography  that  can  also  be  a  potentially  important  modality  for  imaging  of  the  breast 
cancer. 

2.1  Development  of  efficient  linear  reconstruction  algorithms 

In  breast  imaging  applications  of  ultrasound,  the  first-order  Bom  or  Rytov  approximations 
are  typically  not  valid,  and  consequently,  a  linearized  Helmholtz  equation  may  not  accu¬ 
rately  describe  the  relationship  between  the  measured  scattered  wavefield  and  the  scattering 
tissue  [10-18].  However,  many  nonlinear  reconstruction  algorithms  that  account  for  multiple¬ 
scattering  effects,  including  the  ones  we  discuss  below  [19],  involve  the  recursive  application 
of  a  linear  reconstruction  algorithm  (that  assumes  the  validity  of  the  first-order  Bom  or  Rytov 
approximation.)  It  is  therefore  very  important  to  develop  computational  efficient  and  numer¬ 
ically  robust  linear  reconstruction  algorithms  for  DT.  Below,  we  discuss  our  results  on  this 
salient  topic. 

2.1.1  Development  of  3D  efficient  reconstruction  algorithms 

In  2D  DT,  the  filtered  backpropagation  (FBPP)  algorithm  is  [1]  widely  used  for  image  recon¬ 
struction  and  is  generally  regarded  as  being  more  accurate  than  direct  Fourier  reconstruction 
approaches.  However,  objects  such  as  the  female  breast  are  inherently  three-dimensional  and 
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must  be  reconstructed  using  fully  3D  reconstruction  algorithms  in  order  to  avoid  significant 
artifacts  and  a  loss  of  quantitative  accuracy.  We  developed  and  evaluated  novel  reconstruc¬ 
tion  algorithms  for  3D  DT,  referred  to  as  the  esstimation-combination  (E-C)  reconstruction 
algorithms,  that  effectively  solve  the  (fully)  3D  DT  reconstruction  problem  by  performing  a 
series  of  2D  Radon  transform  inversions  [20].  This  greatly  reduces  the  large  computational 
load  that  is  generally  required  by  any  other  3D  DT  reconstruction  technique  such  as  the  3D 
FBPP  algorithm.  This  is  vitally  important  for  the  development  of  computationally  tractable 
3D  nonlinear  algorithms  that  involve  the  recursive  application  of  a  3D  linear  reconstruction 
algorithm.  We  also  demonstrated  that,  in  the  presence  of  data  noise,  there  is  redundant  infor¬ 
mation  contained  in  the  3D  DT  data  function  that  can  be  exploited  by  the  3D  E-C  algorithms 
to  reduce  the  variance  of  the  reconstructed  image  [19]. 

2.1.2  Development  of  minimum-scan  reconstruction  algorithms 

In  many  applications  of  tomographic  imaging  it  is  desirable  to  minimize  the  angular  range 
over  which  the  measurement  data  are  acquired.  This  reduces  the  time  necessary  to  collect 
the  measurement  data,  which  can  reduce  artifacts  due  to  patient  motion.  Furthermore,  in  cer¬ 
tain  situations  it  may  not  be  experimentally  possible  to  collect  data  over  a  complete  2ir  range. 
We  demonstrated  that  a  minimal-scan  data  set  acquired  using  view  angles  only  in  [0,  <f>min ) 
contains  all  the  information  necessary  to  reconstruct  exactly  a  2D  scattering  object  function, 
where  it  <  <  Ztt/2  is  a  function  of  the  measurement  geometry.  Based  on  this  ob¬ 

servation,  we  developed,  investigated,  and  numerically  implemented  minimal-scan  FBPP  and 
E-C  reconstruction  algorithms  for  2D  DT  that  can  exactly  reconstruct  the  scattering  object 
function  from  the  minimal  scan  data  set.  Prior  to  our  work,  all  reconstruction  algorithms 
for  DT  required  a  full  2ir  worth  of  angular  measurements  to  reconstruct  an  accurate  image. 
We  numerically  demonstrated  that  the  minimal-scan  E-C  reconstruction  algorithms  were  less 
susceptible  to  the  effects  of  data  noise  and  inconsistencies  than  were  the  minimal-scan  FBPP 
reconstruction  algorithms.  We  also  generalized  this  work  to  2D  DT  using  the  fan-beam  geom¬ 
etry  and  revealed  a  novel  relationship  between  the  maximum  scanning  angle  and  achievable 
image  resolution.  This  work  may  provide  useful  insights  into  the  development  of  minimum- 
scan  reconstruction  algorithms  for  3D  DT  that  can  be  used  for  breast  imaging  [21]. 

2.2  Development  of  efficient  nonlinear  reconstruction  algorithms 

In  certain  situations  of  the  breast  imaging,  the  first-order  Bom  or  Rytov  approximations  may 
not  be  valid.  Consequently,  a  linearized  Helmholtz  equation  may  not  accurately  describe  the 
relationship  between  the  measured  scattered  wavefield  and  the  scattering  object,  and  nonlinear 
algorithms  are  necessary  for  obtaining  accurate  images.  We  proposed  to  develop  efficient 
nonlinear  reconstruction  algorithms  for  UDT. 


5 


2.2.1  Development  of  2D  efficient  nonlinear  reconstruction  algorithms 

Previously  we  described  our  development  and  investigation  of  E-C  reconstruction  algorithms 
for  linear  DT.  We  have  generalized  these  algorithms  to  the  case  where  a  forward  scattering 
model  includes  multiple-scattering  effects.  Two  forward  scattering  models  were  utilized  that 
captured  higher-order  terms  in  the  Bom  or  Rytov  perturbation  series,  and  are  therefore  po¬ 
tentially  useful  for  modeling  ultrasound  wave  propagation  in  breast  tissue  [22, 23].  For  each 
of  the  two  forward  scattering  models,  we  developed  families  of  nonlinear  E-C  reconstruction 
algorithms  to  solve  the  inverse  problem  [19].  The  nonlinear  E-C  reconstruction  algorithms 
operate  by  relating,  in  2D  Fourier  space  of  the  Radon  transform,  the  n-th  order  perturbation 
of  the  measured  data  function  to  the  n-th  order  perturbation  of  the  scattering  object  function. 
The  algorithms  are  recursive  in  the  sense  that  calculation  of  the  n-th  order  perturbation  of 
the  object  function  requires  knowledge  of  the  (n-l)-th  order  perturbation.  The  computational 
efficiency  of  the  E-C  algorithms  is  therefore  very  relevant  to  this  problem.  We  also  identi¬ 
fied  consistency  conditions  for  the  nonlinear  imaging  models  employed  by  the  two  families  of 
nonlinear  E-C  algorithms.  For  both  imaging  models,  the  consistency  conditions  for  linear  DT 
were  contained  as  special  cases. 

2.2.2  Development  of  3D  efficient  nonlinear  reconstruction  algorithms 

Although  we  have  been  largely  successful  in  the  theoretical  development  of  computationally 
efficient  nonlinear  algorithms  for  2D  UDT,  the  applicability  of  such  algorithms  can  be  re¬ 
strictive  because  the  multi-scattering  effect  in  the  breast  imaging  is  generally  3D  in  nature. 
Therefore,  we  have  also  investigated  3D  nonlinear  reconstruction  algorithms.  Our  strategy  for 
the  development  of  3D  nonlinear  reconstruction  algorithms  is  similar  to  that  for  2D  nonlinear 
reconstruction  discussed  above.  Specifically,  we  proposed  to  investigate  the  two  mentioned 
forward  models  in  3D  and  to  use  the  perturbation  series  for  the  inversion.  The  inversion  of 
the  solution  at  each  perturbative  order  will  be  accomplished  through  the  use  of  our  devel¬ 
oped  linear  3D  E-C  algorithms  for  improving  the  computational  efficiency.  We  are  currently 
continuing  the  development  of  such  3D  perturbative  nonlinear  algorithms. 

2.3  Implementation  and  evaluation  of  the  proposed  algorithms 

We  have  implemented  the  proposed  linear  algorithms  and  nonlinear  algorithms  and  evaluated 
them  by  use  of  computer  simulated  data  and  real  data. 

2.3.1  Implementation  and  evaluation  of  linear  reconstruction  algorithms 

We  have  implemented  the  linear  E-C  reconstruction  algorithms  and  investigated  their  noise 
properties  by  using  a  large  number  of  computer  simulated  data  sets.  Through  our  simula¬ 
tion  studies,  we  have  demonstrated  that  it  is  possible  to  achieve  a  bias-free  reduction  of  the 
statistical  variation  in  the  reconstructed  object  function  by  utilizing  complementary  statistical 
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information  inherent  in  the  scattered  data.  (The  use  of  an  explicit  smoothing  operation  gen¬ 
erally  introduces  bias  in  the  reconstructed  scattering  object  function.)  We  have  quantitatively 
demonstrated  that  the  E-C  algorithms  are  less  susceptible  to  data  noise  and  other  finite  sam¬ 
pling  effects  than  are  the  corresponding  FBPP  algorithms.  This  result  is  consistent  with  the 
observation  that  the  FBPP  algorithms  involve  more  complicated  numerical  operations  (e.g., 
backpropagation)  than  do  the  E-C  algorithms,  which  may  amplify  the  propagation  of  noise 
and  errors  into  the  reconstructed  image.  Using  simulated  strongly  scattering  data,  we  have 
demonstrated  that  the  E-C  algorithms  are  less  susceptible  to  modeling  errors  due  to  viola¬ 
tion  of  first-order  scattering  approximations.  These  same  results  have  been  verified  for  the 
minimum-scan  DT  problem. 

2.3.2  Implementation  and  evaluation  of  nonlinear  reconstruction  algorithms 

Using  simulated  strongly  scattering  data,  we  have  started  to  numerically  investigate  nonlinear 
reconstruction  algorithms  for  2D  DT.  As  described  in  Section  2.2.1,  our  nonlinear  reconstruc¬ 
tion  algorithms  utilize  a  forward  scattering  operator  that  assumes  the  validity  of  a  higher-order 
Bom  or  Rytov  terms.  An  accurate  numerical  implementation  of  the  forward  scattering  oper¬ 
ator  is  critical  for  obtaining  accurate  reconstructions  using  our  algorithms.  In  a  preliminary 
study,  we  have  encountered  difficulty  in  achieving  an  accurate  numerical  implementation  of 
this  operator.  The  forward  scattering  models  employed  by  our  families  of  nonlinear  algorithms 
involve  an  integration  over  a  complex  frequency  variable,  which  is  not  computable  in  practice. 
Accordingly,  numerical  inaccuracies  were  introduced  by  truncating  the  limits  of  integration, 
which  we  observed  to  introduce  a  severe  degradation  in  the  reconstructed  image  quality.  We 
are  currently  investigating  methods  for  mitigating  the  effects  of  the  integration  truncation  used 
by  the  forward  scattering  operator. 

2.4  Development  and  evaluation  of  reconstruction  algorithms  for  reflec¬ 
tivity  tomography 

Reflectivity  tomography  has  been  applied  to  numerous  biomedical  and  non-destructive  test 
imaging  problems  [24-27].  It  has  a  strong  relationship  to  UDT  and  can  be  a  potential  useful 
technique  for  imaging  the  breast  cancer.  The  task  in  reflectivity  tomography  is  to  reconstruct 
from  the  measured  data  a  function  describing  the  reflectivity  distribution  within  the  breast. 
It  has  been  generally  considered  that  accurate  images  can  be  reconstructed  only  from  full 
scan  data  over  In r.  Recently,  we  have  investigated  and  evaluated  image  reconstruction  from 
minimum-scan  data  in  reflectivity  tomography.  Using  the  so-called  potato-peeler  perspective 
that  we  developed,  we  showed  that  accurate  images  can  be  reconstructed  from  minimum-scan 
data  in  reflectivity  tomography.  We  also  performed  quantitative  studies  by  use  of  computer 
simulated  data,  and  the  results  in  such  studies  validated  our  theoretical  results  for  image  re¬ 
construction  in  minimum-scan  reflectivity  tomography. 
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3  Key  research  accomplishments 

•  We  have  developed  and  evaluated  computationally  efficient  3D  linear  reconstruction 
algorithms  that  are  more  than  100  times  faster  than  the  conventional  3D  FBPP  algorithm. 

•  We  have  investigated,  developed,  and  evaluated  algorithms  for  image  reconstruction 
from  minimum-scan  data  in  UDT  with  plane  wave  sources. 

•  We  have  investigated,  developed,  and  evaluated  algorithms  for  image  reconstruction 
from  full-scan  and  minimum-scan  data  in  UDT  with  fan-beam  wave  sources. 

•  We  have  developed  and  evaluated  computationally  efficient  3D  linear  reconstruction 
algorithms  for  UDT  with  spherical  wave  sources. 

•  We  have  developed  computationally  efficient  2D  nonlinear  reconstruction  algorithms 
for  UDT  with  plane  wave  sources. 

•  We  have  investigated  theoretically  the  development  of  computationally  efficient  nonlin¬ 
ear  reconstruction  algorithms  for  3D  UDT. 

•  We  have  developed  computer  codes  that  implement  the  proposed  linear  and  nonlinear 
reconstruction  algorithms. 

•  We  have  evaluated  the  developed  linear  and  nonlinear  reconstruction  algorithms  by  use 
of  computer  simulated  and  experimental  data. 

•  We  have  developed  and  evaluated  reconstruction  algorithms  for  short-scan  reflectivity 
tomography. 

4  Reportable  outcomes 

Peer-Reviewed  Original  Articles 

1.  M.  Anastasio  and  X.  Pan:  Full-  and  minimal-scan  reconstruction  algorithms  for  fan- 
beam  diffraction  tomography,  Appl.  Opt.,  40, 3334-3345, 2001. 

2.  X.  Pan  and  M.  Anastasio:  On  a  limited-view  reconstruction  problem  in  wavefield  to¬ 
mography,  IEEE  Trans.  Med.  Imaging.,  21, 413-416,  2002. 

3.  M.  Anastasio  and  X.  Pan:  Numerically  robust  minimal-scan  reconstruction  algorithms 
for  diffraction  tomography  via  Radon  transform  inversion,  Int.  J.  Imag.  Sys.  Tech.,  (in 
press)  2002. 

4.  M.  Anastasio  and  X.  Pan:  An  improved  reconstruction  algorithm  for  3D  diffraction 
tomography  with  spheric-wave  sources,  IEEE  Trans.  Biomed.  Eng.,  (submitted),  2002. 
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5.  X.  Pan,  Y.  Zou,  and  M.  Anastasio:  Data  redundancy  and  reduced-scan  reconstruction  al¬ 
gorithms  in  reflectivity  tomography,  IEEE  Trans.  Image  Processing,  (submitted),  2002. 

6.  X.  Pan,  Y.  Zou,  and  M.  Anastasio:  Image  reconstruction  of  reflectivity  from  short  scan 
data,  CD  Proc.  of  International  Symposium  on  Biomedical  Imaging,  Washington  D.  C., 
2002. 

Ph.D.  Thesis 

1.  M.  Anastasio:  Development  and  analysis  of  iamge  reconstruction  algorithms  in  diffrac¬ 
tion  tomography,  The  University  of  Chicago,  2001. 

Peer-Reviewed  Proceedings  Articles 

1.  X.  Pan  and  M.  Anastasio:  Minimal-scan  reconstruction  algorithms  for  fan-beam  diffrac¬ 
tion  tomography  and  their  analogy  to  halfscan  fan-beam  CT,  IEEE  Medical  Imaging 
Conference  Record  (CD),  2001. 

2.  M.  Anastasio  and  X.  Pan:  Development  and  evaluation  of  minimal-scan  reconstruction 
algorithms  for  diffraction  tomography,  Proc.  SPIE,  (in  press),  2001. 

Abstracts  and  Presentations 

1.  X.  Pan  and  M.  Anastasio:  Reconstruction  algorithms  in  diffraction  tomography.  Ad¬ 
vanced  Light  Source,  Lawrence  Berkeley  National  Laboratory,  California,  (Host:  Dr. 
Malcolm  Howells),  October  29,  2001. 

2.  M.  Anastasio,  Y.  Zou,  and  X.  Pan:  Reflectivity  tomography  using  temporally  truncated 
data.  The  2nd  Joint  Meeting  of  the  IEEE  Engineering  in  Medicine  and  Biology  Society 
and  Biomedical  Engineering  Society,  (submitted,)  Houston,  2002. 

3.  Y.  Zou,  M.  Anastasio,  and  X.  Pan:  Data  truncation  and  the  exterior  problem  in  reflection¬ 
mode  tomography,  IEEE  Medical  Imaging  Conference,  (submitted,)  Norfolk,  2002. 


5  Conclusion 

Ultrasonic  diffraction  tomography  is  a  potentially  important  technique  for  imaging  of  the 
breast  cancer.  In  this  project,  we  have  investigated,  developed,  and  evaluated  computation¬ 
ally  efficient  and  statistically  optimal  algorithms  for  accurate  reconstruction  of  UDT  images 
that  may  find  applications  in  UDT  imaging  of  breast  cancer.  In  the  last  year,  we  have  made 
contributions  to  UDT  research,  as  summarized  above.  Our  research  on  UDT  have  addressed 
numerous  scientific  and  engineering  problems  involved  in  UDT  image  reconstruction.  These 
results  are  necessary  in  making  UDT  a  viable  medical  imaging  technique  for  imaging  breast 
cancer. 
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Full-  and  minimal-scan  reconstruction  algorithms 
for  fan-beam  diffraction  tomography 


Mark  A.  Anastasio  and  Xiaochuan  Pan 


Diffraction  tomography  (DT)  is  a  tomographic  inversion  technique  that  reconstructs  the  spatially  variant 
refractive-index  distribution  of  a  scattering  object.  In  fan-beam  DT,  the  interrogating  radiation  is  not 
a  plane  wave  but  rather  a  cylindrical  wave  front  emanating  from  a  line  source  located  a  finite  distance 
from  the  scattering  object.  We  reveal  and  examine  the  redundant  information  that  is  inherent  in  the 
fan-beam  DT  data  function.  Such  redundant  information  can  be  exploited  to  reduce  the  reconstructed 
image  variance  or,  alternatively,  to  reduce  the  angular  scanning  requirements  of  the  fan-beam  DT 
experiment.  We  develop  novel  filtered  backpropagation  and  estimate- combination  reconstruction  al¬ 
gorithms  for  full-scan  and  minimal-scan  fan-beam  DT.  The  full-scan  algorithms  utilize  measurements 
taken  over  the  angular  range  0  <  4>  <  2tt,  whereas  the  minimal-scan  reconstruction  algorithms  utilize 
only  measurements  taken  over  the  angular  range  0  s  <j>  <  4>mm>  where  tt  ^  cf)min  <  3tt/2  is  a  specified 
function  that  describes  the  fan-beam  geometry.  We  demonstrate  that  the  full-  and  minimal-scan  fan- 
beam  algorithms  are  mathematically  equivalent.  An  implementation  of  the  algorithms  and  numerical 
results  obtained  with  noiseless  and  noisy  simulated  data  are  presented.  ©  2001  Optical  Society  of 
America 
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1.  Introduction 

Diffraction  tomography  (DT)  is  an  inversion  scheme 
that  can  be  used  for  obtaining  the  spatially  variant 
refractive-index  distribution  of  a  scattering  object. 
Applications  of  DT  can  be  found  in  various  scientific 
fields  such  as  medical  imaging,1'2  nondestructive 
evaluation  of  materials,3*4  and  geophysics.5*6  Unlike 
the  x  rays  used  in  computed  tomography  (CT),  the 
optical  or  acoustical  wave  fields  employed  in  DT  do 
not  generally  travel  along  straight  lines,  thus  pre¬ 
cluding  the  use  of  the  geometrical  optics  approxima¬ 
tion.  Therefore  a  wide  variety  of  techniques  that  are 
suitable  for  reconstruction  of  CT  images  cannot  be 
used  directly  for  reconstruction  of  diffraction  tomo¬ 
graphic  images.  CT  can  be  viewed  as  a  limiting  case 
of  DT,  in  which  the  frequency  of  the  probing  radiation 
tends  toward  infinity. 

A  vast  majority  of  the  algorithm  development  ef- 
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forts  in  DT  have  utilized  the  classic  scanning  geom¬ 
etry,7  which  assumes  that  the  interrogating  radiation 
is  plane  wave  and  the  transmitted  scattered  wave 
field  is  measured  in  a  plane  (or  in  two  dimensions, 
along  a  line)  on  the  opposite  side  of  the  scattering 
object.  This  geometry  is  analogous  to  the  parallel- 
beam  geometry  of  x-ray  CT.  In  many  practical  sit¬ 
uations,  however,  the  interrogating  radiation  may  be 
not  plane  wave  but  rather  produced  by  a  fine  source 
located  a  finite  distance  from  the  scattering  object. 
We  refer  to  this  configuration  as  the  fan-beam  geom¬ 
etry  of  DT,8  which  is  somewhat  analogous  to  the  two- 
dimensional  (2D)  fan-beam  geometry  of  CT. 

The  Born  and  Rytov  approximations9  are  weak- 
scattering  approximations  that  effectively  linearize 
the  inhomogeneous  Helmholtz  and  Ricatti  equations, 
respectively.  The  relative  merits  of  the  Bom  and 
Rytov  approximations  in  the  context  of  DT  have  been 
widely  explored  in  the  literature.10*11  Under  weak- 
scattering  conditions  it  is  customary  and  useful  in  DT 
to  invoke  the  Bom  or  Rytov  approximation  that  per¬ 
mits  the  derivation  of  the  Fourier  diffraction  projec¬ 
tion  (FDP)  theorem.  The  FDP  theorem  relates  the 
one-dimensional  (ID)  Fourier  transform  of  the  mea¬ 
sured  scattered  data  to  the  2D  Fourier  transform  of 
the  scattering  object.  For  2D  DT  employing  plane- 
wave  illumination  and  the  classic  scan  configuration, 
Devaney7  utilized  the  FDP  theorem  to  develop  the 
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well-known  filtered  backpropagation  (FBPP)  algo¬ 
rithm,  which  can  be  viewed  as  a  generalization  of  the 
conventional  filtered  backprojection  (FBPJ)  algo¬ 
rithm  of  x-ray  CT.  Alternative  families  of  plane- 
wave  DT  reconstruction  algorithms,  referred  to  as 
estimate- combination  (E-C)  algorithms  and  general¬ 
ized  FBPP  algorithms,  have  been  developed12  and 
investigated.13'14  The  family  of  plane-wave  E-C  al¬ 
gorithms  effectively  operates  by  transforming  (in  2D 
Fourier  space)  the  DT  problem  into  a  2D  Radon 
transform  problem  that  can  be  efficiently  and  stably 
inverted  by  use  of  the  FBPJ  algorithm.  The  family 
of  plane-wave  generalized  FBPP  algorithms  recon¬ 
structs  the  image  directly  from  the  DT  data  function 
and  includes  the  FBPP  algorithm  as  a  special  mem¬ 
ber.  Both  the  generalized  FBPP  and  E-C  algorithms 
generally  require  scattered  data  measured  from  view 
angles  in  [0,  2tt)  to  perform  an  exact  reconstruction  of 
a  complex-valued  scattering  object.  Accordingly,  we 
refer  to  these  algorithms  as  being  full-scan  algo¬ 
rithms. 

Previously  we  showed15  that,  in  plane-wave  DT 
that  employs  the  2D  classic  scanning  geometry,  a 
minimal-scan  data  set  acquired  by  use  of  view  angles 
only  in  [0,  4>mm  -  3tt/2]  contains  all  the  information 
necessary  for  exact  reconstruction  of  the  scattering 
object  function.  As  the  frequency  of  the  probing  ra¬ 
diation  tends  toward  infinity,  c|>min  ir,  which  re¬ 
flects  the  well-known  fact  that  measurements 
corresponding  to  <$>  G  [0,  tt]  completely  specify  the  2D 
Radon  transform.  [Of  course,  compactly  supported 
objects  are  mathematically  specified  by  a  sinogram 
p(£,  4>0),  where  <$>0  is  contained  in  any  open  set  [0,  e), 
but  if  p(£,  c|>0)  is  not  continuously  sampled  this  obser¬ 
vation  does  not  yield  stable  reconstruction  algo¬ 
rithms.]  We  subsequently  developed  minimal-scan 
FBPP15  and  minimal-scan  E-C16  algorithms  that 
were  capable  of  reconstructing  the  scattering  object 
function  by  use  of  the  minimal-scan  data  set.  Under 
the  conditions  of  continuous  sampling  and  in  the  ab¬ 
sence  of  noise,  we  demonstrated16  that  the  minimal- 
scan  FBPP  and  E-C  reconstruction  algorithms  were 
mathematically  equivalent  to  the  full-scan  FBPP  and 
E-C  reconstruction  algorithms,  respectively. 

Here  we  reveal  and  examine  the  redundant  infor¬ 
mation  that  is  inherent  in  the  fan-beam  DT  data 
function.  Such  redundant  information  can  be  ex¬ 
ploited  to  reduce  the  noise  in  the  reconstructed  image 
or,  alternatively,  to  reduce  the  angular  scanning  re¬ 
quirements  of  the  fan-beam  DT  experiment.  We  de¬ 
velop  novel  E-C  and  FBPP  reconstruction  algorithms 
for  full-scan  and  minimal-scan  fan-beam  DT.  We 
demonstrate  that  the  minimal-scan  algorithms, 
which  utilize  measurements  taken  over  the  angular 
range  0  <  4>  <  <f>min,  where  tt  <  <J>min  ^  3tt/2,  are 
mathematically  equivalent  to  their  full-scan  counter¬ 
parts  that  utilize  measurements  over  the  full  angular 
range  0  <  <|>  <  2tt.  An  implementation  of  the  algo¬ 
rithms  and  numerical  results  obtained  with  noiseless 
and  with  noisy  simulated  data  are  presented. 


/ 

Measured  data  / 


Fig.  1.  Fan-beam  scanning  geometry  of  2D  DT.  The  interrogat¬ 
ing  cylindrical  wave  propagates  along  the  rj  axis,  and  the  scattered 
wave  field  is  measured  along  the  line  q  =  /.  We  obtain  full-scan 
and  minimal-scan  data  sets  by  varying  the  measurement  angle  <j> 
between  0  and  2tt  or  between  0  and  <J>min  [see  Eq.  (29)],  respectively. 
We  assume  that  S  and  D  are  much  larger  than  the  transverse 
dimension  of  the  scattering  object. 

2.  Background 

Here  we  briefly  review  the  geometry  and  approxima¬ 
tions  that  are  traditionally  employed  in  fan-beam  DT, 
as  described  in  the  pioneering  work  of  Devaney.8  A 
table  of  the  acronyms  used  in  this  manuscript  is  in¬ 
cluded  in  Appendix  A. 

A.  Fan-Beam  Diffraction  Tomography 

The  classic  scanning  configuration  of  fan-beam  DT  is 
shown  in  Fig.  1.  The  fixed  coordinate  system  (x,y), 
the  rotated  coordinate  system  (£,  -q),  and  the  usual 
polar  coordinates  (r,  0)  are  related  by  x  =  r  cos  0,  y  = 
r  sin  0,  ^  =  x  cos  $  +  y  sin  cj)  =  r  cos  (4>  -  0),  and  r\  = 
— x  sin  c|>  +  y  cos  (J >  =  -r  sin(<|>  -  0).  The  scattering 
object  is  illuminated  by  a  monochromatic  cylindrical- 
wave  source  located  at  the  position  tj  =  -S  on  the  rj 
axis,  emitting  a  wave  field  of  the  form 

/,  n  „  exp(j2iTVo|r  -  Sp|) 

40  =  u° - iT^sli - 

exp{j2irv0[f  +  QS  +  Z))2]1/2} 

0  [£2  +  (S  +  D)2]1'2  ’ 

where  U0  is  the  complex  amplitude,  k  =  2ttv0  is  the 
wave  number,  p  is  a  unit  vector  pointing  along  the 
positive  T|  axis,  and  D  is  the  distance  of  the  detector 
surface  from  the  center  of  rotation.  The  incident 
wave  field  u ,(£,  <(>)  could  represent  a  pressure  field  in 
acoustical  applications  or  a  scalar  electromagnetic 
field  in  optical  applications ,  for  example.  From  mea¬ 
surements  of  the  scattered  wave  field  obtained  along 
the  £  axis  oriented  at  a  measurement  angle  $  at  a 
distance  r\  =  D  from  the  origin,  one  seeks  to  recon¬ 
struct  the  scattering  object  function  a(r).  The  un¬ 
derlying  physical  property  of  the  scattering  object 
that  is  mapped  in  DT  is  the  refractive-index  distri- 
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bution  n( r),  which  is  related  to  the  scattering  object 
function  by  the  equation  a(r)  =  n2(r)  -  1. 

Let  «(£,  <}>)  denote  the  measured  total  wave  field 
along  the  line  r)  =  D,  as  shown  in  Fig.  1.  The  scat¬ 
tered  data  are  given  by 

uM,  <i>)  =  «(€,  <|>)  -  u,(£,  <\>),  (2) 

which  can  be  treated  as  a  measurable  data  function 
because  u(£,  <t>)  and  <j>)  can  be  measured.  There¬ 
fore  we  can  introduce  a  modified  data  function  M(vm, 
<J>),  which  is  given  by 


M(vm,  <|>)  = - j  v'  exp[-j'2Tt(v'  -  v0)D] 

HVo 


uM,  <1 >) 

Ui(£,  4>). 


(3) 


Fig.  2.  The  fan-beam  FDP  theorem  states  that  M(vm ,  4>)  is  equal 
to  the  2D  Fourier  transform  of  a(r)  along  the  semielliptic  arc  AOB 
that  has  semiaxes  equal  to  v0  and  vQ/\. 


where 


=  (l/2ir) 


fi(€)exp(-y'2Trvm0d^ 


B.  Fan-Beam  Fourier  Diffraction  Projection  Theorem 

In  plane-wave  DT,  the  FDP  theorem17  relates  the 
modified  data  function  to  the  2D  Fourier  transform  of 
the  scattering  object  and  can  be  viewed  as  a  general¬ 
ization  of  the  Fourier  slice  theorem  of  conventional 
x-ray  CT.  The  FDP  theorem  is  valid  under  condi¬ 
tions  of  weak  scattering  and  plane-wave  illumina¬ 
tion.  To  establish  the  FDP  theorem  for  the  fan- 
beam  DT  geometry  it  is  necessary  to  assume  the 
weak-scattering  approximation  and  the  so-called 
paraxial  approximation,8  which  is  a  well-known  ap¬ 
proximation  in  the  optics  literature.  The  paraxial 
approximation  requires  that  both  S  and  D  be  much 
larger  than  the  dimension  size  of  the  scattering  ob¬ 
ject.  This  amounts  to  requiring  that  the  largest  an¬ 
gle  subtended  by  the  object  when  the  object  is  viewed 
from  either  the  source  or  the  measurement  location 
be  much  smaller  than  a  radian. 

Under  the  Bom  and  paraxial  approximations,  Dev- 
aney  derived  the  fan-beam  FDP  theorem,8  which  is 
given  by 


As  displayed  in  Fig.  2,  Eq.  (6)  states  that  the  modified 
data  function,  M(vm,  <[>)>  is  equal  to  a  semielliptical 
slice,  oriented  at  angle  cj>,  through  the  2D  Fourier 
transform  of  the  object  function  a(r).  One  can  also 
derive  the  FDP  theorem  by  employing  the  Rytov  ap¬ 
proximation  instead  of  the  Bom  approximation.  In 
this  case,  Eq.  (6)  remains  unchanged  and  only  Eq.  (3) 
needs  to  be  appropriately  redefined.8 

3.  Full-Scan  Reconstruction  Algorithms  for  Fan-Beam 
Diffraction  Tomography 

First,  we  present  families  of  full-scan  FBPP  and  E-C 
reconstruction  algorithms  for  fan-beam  DT.  These 
fan-beam  algorithms  are  novel  and  contain  the  pre¬ 
viously  developed  families  of  plane-wave  FBPP  and 
E-C  algorithms  as  limiting  cases.  They  will  be  gen¬ 
eralized  to  the  minimal-scan  situation  in  Section  4 
below. 

A.  Fan-Beam  Full-Scan  Estimate-Combination 
Algorithms 

The  Radon  transform18  of  the  scattering  object  func¬ 
tion  a(r,  0)  is  defined  as 


a(r,  0)8[|  -  r  cos(<}>0  -  0)]dr, 

(7) 

where  cf>0  is  the  projection  angle,  £  =  r  cos(4>0  —  0),  and 
tj  =  t  sin(4>0  -  0).  The  2D  Fourier  transform  of p( 
<|>o)  is  defined  as  [strictly  speaking,  P*(vG)  is  the  com- 


M(vm,  <(>)  =  J  |  a(r)exp|  -;2tt|^  £  -  (v'  -  v0)t]  j  jdr  |vm|  <  xv0, 
=  0  |vm|  >  xv0- 
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bination  of  the  ID  Fourier  transform  with  respect  to 
va  and  a  ID  Fourier  series  expansion  with  respect  to 
4>0  of  the  Radon  transform] 


pit  4>o) 


X  exp[-y2'irva£  -jk<$> 0]d£d<j>0, 


(8) 


where  va  is  the  spatial  frequency  with  respect  to  £  and 
the  integer  k  is  the  angular  frequency  index  with 
respect  to  <t>0.  It  is  well  known  that  a(r,  0)  can 
readily  be  reconstructed  exactly  from  its  Radon 
transform  p(£,  <|>0)  [or,  equivalently,  PA(va)]  by  use  of 
a  wide  variety  of  computationally  efficient  and  nu¬ 
merically  stable  reconstruction  algorithms  such  as 
the  FBPJ  algorithm.  Therefore  the  task  of  image 
reconstruction  in  fan-beam  DT  is  tantamount  to  the 
task  of  estimating  Pk(va)-  Furthermore,  because 
Pk(va)  for  va  >:  0  contains  full  knowledge  of  the  Radon 
transform,  one  needs  to  estimate  P*,(va)  from  the  mea¬ 
sured  scattered  data  only  for  va  >  0. 


Comparison  of  Eqs.  (9)  and  (11)  yields  that,  for  |vm| 
^  Xv0> 

Pk(va)  =  y(vm/x2)kMk(vm),  (13) 

provided  that 


From  Eq.  (14)  we  see  that  va  is  real  (that  is,  0  ^  va  ^ 
V0  Vl  +  the)  only  for  |vj  <  xvo- 
In  the  absence  of  data  noise  or  inconsistencies,  one 
can  use  Eq.  (13)  to  obtain  Pk(va)  exactly  from  Mk(vm)f 
which  can  readily  be  obtained  from  the  modified  data 
function.  In  the  presence  of  data  noise  or  inconsis¬ 
tencies,  one  can  use  Eq.  (13)  to  obtain  an  estimate  of 
P*(vB).  For  any  given  0  <  va  <  v0Vl  +  1/x2  »  we 
show  in  Appendix  B  that  four  different  roots  vmi,  i  =  1, 
2,  3,  4,  satisfy  Eq.  (14).  However,  only  two  of  these 
four  roots  correspond  to  real-valued  frequencies,  which 
are  given  by 


vmi  =  ”Vm2  =  vm  =  va  1  - 


2  \  1/2 


4v02/  \V2{1  -  (1  -  x2)(v„2/ 2V02)  +  [1  -  x2(l  -  X2)(va2/v02)]1/2}. 


1/2 


(15) 


From  Eqs.  (7)  and  (8)  it  can  be  shown19  that 

/'2tt  /*°° 

Pk(va)  =  (~j)k  «(r) 

Je=0  * r=0 

X  exp(-7^0)Jjfe(2'irvar)rdrd0,  (9) 

where  Jk  indicates  the  ^th- order  Bessel  function  of 
the  first  kind.  Because  M(vmi  <|))  is  a  periodic  func¬ 
tion  of  4>,  it  can  be  expanded  into  a  Fourier  series  with 
expansion  coefficients  given  by 


Mk(vm)  = 


=  -f 

2”J„ 


M(vm,  <J))exp(-j&4>)d()>.  (10) 


Substituting  Eq.  (6)  into  Eq.  (10),  noting  that  f  =  r 
cos(c()  —  0)  and  iq  =  — r  sin(<})  —  0),  and  carrying  out  the 
integration  over  angle  4>  (Ref.  19)  yield 


In  the  plane-wave  case  there  are  only  two  roots,12 
which  one  can  obtain  from  Eq.  (15)  by  letting  x  =  1. 

Therefore,  for  a  given  0  <  v„  <  v0(l  +  l/x2)1/2,  one 
can  obtain  two  estimates  of  PA(va)  from  knowledge  of 
Mk(vm)  at  vml  and  vm2,  namely, 


PM  = 


Mk(vm\)  = 


Mk(vm), 

(16) 


Pk(va)  = 


h 

Mk(vm2) 


-k 


(17) 


In  establishing  Eq.  (17)  we  used  the  property  -y(-vm') 
=  1.  It  can  readily  be  shown  that,  as  the 


02ir  _ _ 

a(r)exp(-jkQ)J k(2nr  JvJ2  -  v/)rdrd0 

.  )=0 

=  0 


\Vm\  -  X^O, 

|vm|  >  Xv0) 


(ID 


where  vm'  =  vm/x2,  =  j(v'  -  v0),  and 


y(vm) 


Vvm'2  ~  Vp.2 

Vm'  +  Vp. 


incident  wave  becomes  planar  (i.e.,  as  x  -*  1),  the 
above  results  become  the  results  in  Ref.  12  for  the 
plane-wave  case. 

In  the  absence  of  noise,  Eqs.  (16)  and  (17)  yield 
(12)  identical  (and  exact)  values  of  Pk(va).  In  the  pres¬ 
ence  of  data  noise  (or  other  inconsistencies),  the  two 
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estimates  of  Pk(va)  are  distinct,  suggesting  that  it 
may  be  beneficial  to  combine  the  two  estimates  lin¬ 
early12  to  form  a  final  estimate  of  PA(va)  as 

Pk(va)  =  0>k(vm)hkMk(vm)] 

+  [1  -  «*(vm)][(-l)S-*M*(-vJ],  (18) 

where  (oA(vm)  is  a  combination  coefficient  that  dic¬ 
tates  how  Eqs.  (16)  and  (17)  are  combined.  This 
strategy  of  linear  combination  has  been  demon¬ 
strated  to  be  useful  in  reducing  the  noise  in  the  re¬ 
constructed  image.12-14  Because  each  selection  of 
t ok(vm)  gives  rise  to  a  particular  final  estimate  P*(va), 
Eq.  (18)  can  be  interpreted  as  an  estimation  method 
for  obtaining  P*(va)  (or,  equivalently,  the  Radon 
transform).  Because  o)A(vm)  may  be  any  complex¬ 
valued  function  of  vm  and  k ,  Eq.  (18),  in  effect,  pro¬ 
vides  infinite  families  of  estimation  methods.  From 
the  estimate  P^(va)  (i.e.,  the  Radon  transform),  one 
can  subsequently  reconstruct  the  image  a(r)  by  use  of 
the  FBPJ  algorithm.  For  simplicity,  the  use  of  Eq. 
(18)  to  estimate  Pk(va)  coupled  with  the  2D  FBPJ 
algorithm  to  reconstruct  a(r)  is  referred  to  as  a  fan- 
beam  full-scan  E-C  reconstruction  algorithm.  As  S 
oo,  we  observe  that  x  1  and  that  the  fan-beam 
full-scan  E-C  reconstruction  algorithms  reduce  to  the 
plane-wave  full-scan  E-C  reconstruction  algorithms 
developed  previously.12 

B.  Fan-Beam  Full-Scan  Filtered  Backpropagation 
Algorithms 

The  fan-beam  full-scan  E-C  algorithms  discussed 
above  first  estimate  Pk(va)  (i.e.,  the  Radon  transform) 
from  the  modified  data  function  Mk(vm)  and  subse¬ 
quently  reconstruct  the  image  by  inverting  the  esti¬ 
mated  Radon  transform.  Below,  we  develop 
algorithms  that  reconstruct  images  directly  from  the 
modified  data  function.  Using  Eq.  (9),  one  can  di¬ 
rectly  express  the  object  function  in  terms  of  the  es¬ 
timate  Pk{va)  as 

o(r)=.2ir  2/  f"  Pk{va) 

k~  ”  *  Vo=0 

X  exp(jkQ)Jk(2TTvar)vadva.  (19) 


Substituting  this  result  into  Eq.  (20),  and  noting  that 
»■'  =  Vm/X2,  yields 


00  /*XV0  \f 

a(r)  =  it  2  /  exp(j^e)  —  [(1  -  xV 

*—  Jvm=0  X  V 

+  X2v0]WvJ7*M*(vJ  +  [1  -  w*(vm)] 

X  ( - 1)*7  ~kMk{  - vm) } 

X  Jk[2i:(vm’2  -  v^r)1/2]dvm.  (22) 

To  reduce  Eq.  (22)  to  the  form  of  the  FBPP  algo¬ 
rithm,  we  assume  that  uk(vm)  +  (x>k(~ vm)  =  1.  Using 
the  integral  identities12 

(±)kjky(vm')±k  exp(jkQ)Jk[2n(vm'2  -  v1>)1/2] 


»2it 

% 

JH  ±j 2lTVm  -j  +  2lTV  T| 

exp 

Jo 

X  J 

and  defining 


00 


4*)  =  2  exp(j^4>)wA(vm)M*(vm) 
^=-00 


yield  for  Eq.  (22) 


(24) 


1  f2™  /Vo  Iv  I 

«(•*)  =  2  ~T7  [(!  -  xV'  +  X2Vo]M<“)(vm,  <|>) 

Jo 


X  exp 


2  +  2™^ 

X 


d<j)dvm. 


(25) 


Equation  (25)  describes  a  family  of  fan-beam  full-scan 
FBPP  algorithms  [indexed  by  the  choice  of  a>^(vm)], 
which  becomes  the  family  of  FBPP  algorithms  for 
plane-wave  full-scan  DT12  when  x  1.  In  particu¬ 
lar,  when  wA(vm)  —  1/2,  Eq.  (25)  corresponds  to  the 
fan-beam  FBPP  algorithm  suggested  by  Devaney.8 
Because  the  derivation  of  the  family  of  fan-beam  full- 
scan  FBPP  algorithms  was  based  on  the  family  of 
fan-beam  full-scan  E-C  algorithms  [Eq.  (20)],  the  E-C 
and  FBPP  algorithms  are  mathematically  equiva¬ 
lent.  However,  as  will  be  demonstrated  below,  these 
FBPP  and  E-C  algorithms  respond  differently  to  data 
noise  and  other  experimental  errors. 


Substituting  Eq.  (18)  into  Eq.  (19)  yields 

"  f  v0  y/l  +  l/x2 

“(»*)  =  2-rr  2  J  exp (jkQ)  Wk(vm)y  Mk{vm) 

Jv0— o 

+  [1  -  ^(vJ](-l)*7-*M*(-vmM(27rv0r)v0dv0. 

(20) 

Using  the  relationship  between  va  and  vm  in  Eq.  (14), 
one  can  show  that,  for  0  <  va  <  v0Vl  +  1/x2, 

=  Vm  [(1  -  xV'  +  (21) 


4.  Minimal-Scan  Reconstruction  Algorithms  for 
Fan-Beam  Diffraction  Tomography 

As  discussed  above,  the  redundant  information  con¬ 
tained  in  the  scattered  data  can  be  employed  for  re¬ 
ducing  the  image  noise.  Such  information  can  also 
be  used  for  reducing  the  angular  scanning  require¬ 
ments  in  a  DT  experiment. 

A.  Consistency  Conditions  and  the  Fan-Beam  Data 
Space 

According  to  the  fan-beam  FDP  theorem  in  Eq.  (6), 
the  modified  data  function  M(vm,  4>)  satisfies  the  con¬ 
sistency  condition 

M(vm,  <|>)  =  4>  +  TT  -  2a),  (26) 
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Fig.  3.  Complete  data  space  tW  =  .s4US&U<:€U2>  contains  data 
from  the  view  angles  in  [0,  2t it].  Subset 

obtained  from  the  view  angles  in  [0,  <f>min]  is  called  the  minimal- 
complete  data  set.  The  minimal-complete  data  set  contains  all 
the  information  necessary  for  exact  reconstruction  of  the  scattering 
object  function. 


where 


sin  a  =  sgn(vm) 


(v'  ~  vo)2 

(vm2/x2)  +  (v'  -  vo)2. 


1/2 

.  (27) 


Let  <W  =  [|vm|  s  xvo>  0  <  4>  s  2ir]  denote  the  complete 
(or  full-scan)  data  set.  As  shown  in  Fig.  3,  °W  can  be 
divided  into  the  four  subspaces,  si,  95,  c€,  and  9), 
where 


&  =  [|vm|  Xv0,  ,0  <  4>  <  2ot  +  28], 

28  =  [|vm|  ^  xvoj  2a  +  28  <  <)>  <  -tt  4-  2a], 
^  =  [|vj  ^  Xvo>  it  +  2a  <  d>  <  4>min], 

=  [|vm|  ^  Xv0>  4>n.in  ^  <M  2ir], 

where 


sin  8  =  (1  +  l/x2)1/2  ’  tf>min  =  7r  +  28'  (28) 

Using  Eqs.  (26)  and  (27),  one  can  verify  that  infor¬ 
mation  in  subspace  %  is  identical  to  that  in  subspace 
si  and  that  information  in  subspace  S6  is  identical  to 
that  in  subspace  2).  As  shown  in  Fig.  3,  because  the 
boundaries  between  the  subspaces  are  generally 
functions  of  vm  and  <j>  and  because  each  horizontal 
line  in  corresponds  to  a  measurement  acquired  at 
a  particular  view  angle,  the  information  in  subspace 
Sft  cannot  in  practice  be  determined  independently  of 
that  in  subspace  %  and  vice  versa.  We  therefore 
refer  to  the  union,  it  =  ^U2SUc€  =  [|vm|  ^  xvo>  0  — 
4>  —  4>min3>  as  minimal-complete  data  set.  A  plot 
of  <|)min  versus  x  is  shown  in  Fig.  4.  As  x  — >  1, 4>min  — » 
3tt/2,  and  M  reduces  to  the  minimal-complete  data 
set  proposed  previously  for  plane-wave  DT.15  How¬ 
ever,  for  x  <  1,  <t>rain  <  3tt/2,  which  indicates  that  the 


Fig.  4.  Plot  of  <|)min  versus  x  reveals  that,  in  fan-beam  DT,  the 
required  angular  scanning  range  is  less  than  270°. 


angular  scanning  requirements  of  fan-beam  DT  are 
less  restrictive  than  for  plane-wave  DT. 

Figure  5  clearly  demonstrates  that  the  minimal- 
complete  data  set  contains  all  the  information  re¬ 
quired  for  an  exact  reconstruction.  According  to  the 
FDP  theorem,  the  segment  AOB  corresponds  to  a 
semielliptical  slice  through  the  2D  Fourier  transform 
of  a(r),  which  can  be  obtained  from  the  modified  data 
function  M(vm,  c}>).  The  Fourier  space  coverages  pro¬ 
duced  by  the  segments  OA  and  OB  as  4>  varies  from 
0  to  4>min  are  shown  in  Figs.  5(a)  and  5(b),  respec- 


(a)  (b) 


(c) 

Fig.  5.  As  <)>  varies  from  0  to  4>min,  the  two  segments  OA  and  OB 
(Fig.  2)  yield  two  incomplete  coverages  of  the  2D  Fourier  space  of 
the  object  function,  which  are  shown  in  (a)  and  (b),  respectively. 
Superimposing  the  two  incomplete  coverages  in  (a)  and  (b),  one 
obtains  a  complete  coverage,  as  shown  in  (c),  of  the  2D  Fourier 
space  of  the  object  function. 
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tively.  It  can  be  observed  that  each  of  these  two 
coverages  alone  provides  only  an  incomplete  coverage 
of  the  2D  Fourier  space  of  a(r).  However,  one  can 
superimpose  these  two  incomplete  coverages  in  Figs. 
5(a)  and  5(b)  to  obtain  complete  coverage  of  the  2D 
Fourier  space  of  a(r),  as  shown  in  Fig.  5(c). 

The  redundant  information  contained  in  subspaces 
si  and  ^  of  the  minimal-complete  data  set  needs  to  be 
normalized  properly  before  or  during  the  reconstruc¬ 
tion  procedure.  Let  40  denote  the  minimal 

scan  data,  where  M(m)(vm,  4)  =  M(vm,  4>)  for  0  <  4  < 
4>min  and  M<m)(vm,  <J>)  =  0  for  <{)min  <  cj>  <  2tt.  Con- 
sider  a  weighted  data  set  M'(vm,  4),  defined  as 

4>)  =  w(vm,  4>)M<m\vm,  4),  (29) 

where  w(vm ,  4)  can  be  a  function  of  vm  and  4>  that 
satisfies 

w(vm,  <|>)  +  w(-vm,  c|>  +  -it  -  2a)  =  1  (30a) 

in  complete  data  space  °W, 

w(vm,  4)  =  1  (30b) 

in  subspace  SS,  and 

w(vm,  4)  =  0  (30c) 

in  subspace  2L  One  can  choose  different  w(vm ,  4)  in 
subspaces  si  and  ^  as  long  as  these  w(ym,  4)  satisfy 
Eq.  30(a).  In  the  numerical  examples  that  follow, 
we  used  the  explicit  form  for  w(vm,  4)  given  by  Eq. 
(38).  We  can  now  readily  obtain  minimal-scan  re¬ 
construction  algorithms  for  fan-beam  DT. 

B.  Fan-Beam  Minimal-Scan  Estimate-Combination 
Algorithms  1 

Because  of  Eq.  30(c),  Eq.  (29)  can  also  be  rewritten  as 
M'(vm>  <J»)  =  w(ym,  <|>)M(vm,  c|>).  (31) 

Using  Eqs.  (26)  and  (30a),  one  can  verify  that 

M(vm,  <|>)  =  M'(vm,  <(>)  +  <(>  +  u  -  2a). 

(32) 

Using  Eq.  (32)  in  Eq.  (10),  one  obtains 

Mkiy  J  =  [Mk'(v  J  +  (-l)V^'(-vJ],  (33) 

where 

P2ir 

Mk'(vm)  =  (l/2ir)  exp {-jk$)M'{vm,  <t>)d<4>. 

Jo 

Multiplying  both  sides  of  Eq.  (33)  by  yk  and  noting 
Eq.  (16),  we  conclude  that,  for  \vm\  <  xvo> 

Pk(va)  =  LykMk'(vm)  +  (-1  )ky-kMk’(-vm)l  (34) 

where  va  and  vm  are  related  by  Eq.  (14).  A  fan-beam 
minimal-scan  E-C  algorithm  is  formed  by  use  of  Eq. 
(34)  to  estimate  the  Radon  transform  from  the 
minimal-complete  data  set  acquired  at  measurement 
angles  0  <  4  ~  4>min  and  the  FBPJ  algorithm  to 
reconstruct  the  final  image.  One  can  form  different 


Fig.  6.  The  numerical  phantom  used  in  the  simulation  studies 
comprised  two  concentric  ellipses.  The  fan-beam  FDP  theorem 
was  used  to  calculate  analytically  the  simulated  scattered  field 
data  from  the  phantom. 


fan-beam  minimal-scan  E-C  algorithms  by  specifying 
different  choices  for  w(vm ,  4)  in  Eq.  (31)  that  satisfy 
Eqs.  (30). 

C.  Fan-Beam  Minimal-Scan  Filtered  Backpropagation 
Algorithms 

Using  Eq.  (34)  in  Eq.  (19)  and  using  the  strategy 
outlined  in  Subsection  3.B,  we  can  also  develop  fan- 
beam  minimal-scan  FBPP  algorithms  that  can  recon¬ 
struct  the  object  function  directly  from  the  weighted 
minimal-complete  data  set  that  is  given  by 


nxv°  |VJ  0 

^[(i-xV  +  xV|M'(vm)<f>) 

-YVn  X  V 


X  exp  j*27Tvm  2  +  2tcv  y\ 
'  L  X 


dc|>dv„ 


(35) 


where  4>min  is  a  function  of  x  as  stated  in  Eq.  (28). 
One  can  form  different  fan-beam  minimal-scan  FBPP 
algorithms  by  specifying  different  choices  for  w(vm ,  4>) 
in  Eq.  (31)  that  satisfy  Eqs.  (30).  When  x  =  1,  the 
fan-beam  minimal-scan  FBPP  algorithms  reduce  to 
the  plane-wave  minimal-scan  FBPP  algorithms  de¬ 
rived  previously.15 


5.  Numerical  Results 

We  numerically  investigated  the  fan-beam  full-  and 
minimal-scan  reconstruction  algorithms,  using  sim¬ 
ulated  noiseless  and  noisy  data. 

A.  Data  and  Noise  Model 

We  employed  the  numerical  phantom  composed  of 
two  concentric  ellipses,  as  displayed  in  Fig.  6.  The 
values  of  the  scattering  object  function  that  corre¬ 
spond  to  the  outer  and  inner  ellipses  are  0.0005  and 
0.0001,  respectively.  We  chose  a  fan-beam  geometry 
specified  by  x  “  0.8,  but  our  observations  below  hold 
for  arbitrary  x*  The  FDP  theorem  was  employed  to 
calculate  analytically  the  modified  data  function 
M(vm ,  4)  that  by  means  of  Eq.  (3)  determined  the 
scattered  field  data  us(£,  4).  Therefore  our  simula¬ 
tions  were  designed  to  demonstrate  the  performance 
of  the  reconstruction  algorithms  under  the  condition 
that  the  Bom  and  paraxial  approximations  are  valid. 
The  evaluation  of  the  performance  of  the  algorithms 
when  the  Born  and  paraxial  approximations  are  not 
valid10'11-20  remains  a  topic  for  future  study.  The 
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discrete  complete  data  set  comprised  128  equally 
spaced  measurement  angles  in  [0,  2rr).  The  discrete 
minimal-complete  data  set  comprised  92  equally 
spaced  measurement  angles  in  [0,  4.49]  (or,  equiva¬ 
lently,  [0,  257.3°]).  In  this  way,  both  data  sets  had 
the  same  angular  sampling  increment,  Ac}>  =  2tt/ 
128  «  4,49/92.  The  data  function  M(vm,  <}))  con¬ 
tained  129  evenly  spaced  samples  in  [“Xvo>  Xvol- 
To  simulate  the  effects  of  data  noise,  we  treated  the 
scattered  data  us(£,  <}))  as  a  complex  stochastic  pro¬ 
cess  with  a  real  and  an  imaginary  component,  de¬ 
noted  us(r)  (4,  <\>)  and  4>),  respectively.  (Here, 

boldface  type  for  u  and  a  denotes  a  random  variable.) 
Let  us(r)  =  us(r)  +  Aus(r)  and  us(l)  =  u}l)  +  Aus(l), 
where  us(r)  and  us(l)  are  the  means  of  us(r)  and  us(i). 
respectively.  The  statistics  of  the  deviates  Aus(r' 
and  Aus(l)  are  described  by  the  circular  Gaussian 
model, 


p(Ausw,  Au®)  = 


2ivcri(T 
X  exp 


(A  u*  ^ 

Au®2\l 

2 

+ 

^  jj 

(36) 


where  crr2  and  u?  are  the  variances  of  Aus(r)  (£,  <}))  and 
Aus(l)(£,  4>),  respectively. 

To  study  the  noise  properties  of  the  reconstructed 
images  quantitatively,  we  generated  N  -  250  noisy 
complete  and  minimal-complete  data  sets  by  using 
the  noise  model  in  Eq.  (36)  with  oy  =  a*  =  0.05.  We 
used  the  fan-beam  full-scan  and  minimal-scan  E-C 
and  FBPP  algorithms  to  reconstruct  sets  of  250  noisy 
images  from  these  noisy  data  sets.  The  matrix  size 
of  the  reconstructed  images  was  128  X  128  pixels, 
and  the  wavelength  of  the  incident  radiation  was 
equal  to  2  pixels.  The  local  image  variance  was  cal¬ 
culated  empirically  from  the  N  sets  of  reconstructed 
images  as 


var[a(r)]  = 


N-  1 


I**-* 


N 

2  ai(r) 


12 


(37) 

where  a£(r)  is  the  ith  image  obtained  by  use  of  the 
reconstruction  algorithm  under  investigation. 


of  Mk(vm  =  VV  +  V1)  and  Mk(vm  =  -VvJ  +  v^) 
from  the  sampled  values  of  Mk(vm)  and  Mk(-vm), 
which  we  subsequently  used  in  Eq.  (18)  to  evaluate 
Pk(va).  For  each  value  of  k ,  zero-padding  interpola¬ 
tion  was  employed  to  increase  the  sampling  density 
along  the  vm  axis  of  Mk(vm)  by  a  factor  of  3  to  increase 
the  accuracy  of  the  interpolation  operation.  We  em¬ 
ployed  the  consistency  condition  PA(va )  = 
(-  l)kPk(-va)  to  obtain  the  64  samples  of  Pk(va)  for  va 
spanning  the  range  -v0Vl  +  1/x^  ^  va  <  0.  In  our 
implementation  of  the  FBPJ  algorithm,  an  unapo- 
dized  ramp  filter  was  used.  The  interpolation  nec¬ 
essary  for  aligning  the  backprojected  data  onto  a 
128  X  128  pixel  discrete  image  matrix  was  performed 
by  bilinear  interpolation.  When  Mk(vm)  is  replaced 
by  Ma'(vw),  and  Eq.  (34)  is  employed  in  place  of  Eq. 
(18),  the  above  paragraph  also  describes  our  imple¬ 
mentation  of  the  fan-beam  minimal-scan  E-C  algo¬ 
rithm. 


2.  FBPP  Algorithms 

In  the  full-scan  and  minimal-scan  FBPP  algorithms, 
at  each  measurement  angle  c|>,  M(vm,  <}>)  orM'(vm,  <j>), 
respectively,  was  multiplied  bv  the  depth-dependent 
filter  (|vj/x4v')[(l  -  x2)v'  +  X^olexp^irv^ti]  for  each 
of  128  discrete  values  of  y\.  For  each  value  of  v\,  the 
filtered  data  were  zero  padded  to  ensure  that  the 
pixel  size  of  the  reconstructed  image  matched  the 
pixel  size  of  the  images  reconstructed  by  use  of  the 
full-  and  minimal-scan  E-C  algorithms.  The  inter¬ 
polation  necessary  for  aligning  the  backpropagated 
data  onto  a  128  X  128  pixel  discrete  image  matrix 
was  performed  by  bilinear  interpolation. 

The  fan-beam  minimal-scan  E-C  and  minimal-scan 
FBPP  algorithms  utilized  the  weight  function  w{vm, 
4>)  [see  Eq.  (29)]  given  by 


ui{vm,  <t>) 


«iv.2 

'll  (|> 

sin 

4  (tt/4)  —  a 

i 

[ IT  (3-77/2)  - 

*1 

sin 

[4  (Tr/4)  +  a  J 

10 


0  <  <f  <  25  +  2a 
2S  +  2a  <  cj)  <  tt  +  2a 
7T  +  2a  <  4>  <  <}>min 
<!>min  ^  4>  =£  2ir 

(38) 


B.  Implementation  Details 
L  E-C  Algorithms 

From  the  uniformly  sampled  values  of  the  scattered 
field  ms(£,  <J>),  Mk(vm)  can  be  determined  at  uniformly 
spaced  values  of  vm.  However,  because  of  the  non¬ 
linear  relationship  [Eq.  (14)]  between  va  and  vm,  the 
uniformly  spaced  values  of  vm  at  which  Mk(vm)  is 
sampled  do  not  generally  correspond  to  the  uniformly 
spaced  values  of  va  at  which  one  needs  to  evaluate 
Pk(va)  (to  utilize  the  fast  Fourier  transform  in  the 
FBPJ  algorithm).  For  each  of  65  evenly  spaced  val- 
ues  of  va  spanning  the  range  0  <  va  <  v0Vl  +  l/x*, 
we  used  linear  interpolation  to  determine  the  values 


C.  Results 

From  the  simulated  noiseless  complete  and  minimal- 
complete  data  sets  we  reconstructed  the  phantom, 
using  the  full-  and  minimal-scan  fan-beam  recon¬ 
struction  algorithms.  Figures  7(a)  and  7(b)  show 
the  images  obtained  from  the  full-scan  and  the 
minimal-scan  E-C  reconstruction  algorithms,  respec¬ 
tively.  The  full-scan  E-C  algorithm  was  specified  by 
o)fcC vm)  =  1/2  in  Eq.  (18).  The  images  appear  iden¬ 
tical,  as  is  consistent  with  our  assertion  that  the  full- 
and  minimal-scan  E-C  algorithms  are  mathemati¬ 
cally  equivalent  in  the  absence  of  noise  or  other  er¬ 
rors.  Figures  7(c)  and  7(d)  show  the  images 
obtained  by  use  of  the  full-scan  and  the  minimal-scan 
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(C)  (d) 

Fig.  7.  Images  reconstructed  from  noiseless  data  by  use  of  (a) 
full-scan  E-C,  (b)  minimal-scan  E-C,  (c)  full-scan  FBPP,  and  (d) 
minimal-scan  FBPP  reconstruction  algorithms. 

FBPP  reconstruction  algorithms,  respectively. 
Again,  the  images  appear  identical,  consistent  with 
our  assertion  that  the  full-  and  minimal-scan  FBPP 
algorithms  are  mathematically  equivalent  in  the  ab¬ 
sence  of  noise  or  other  errors.  As  expected,  it  is  also 
observed  that  the  images  reconstructed  with  the  full- 
and  minimal-scan  E-C  algorithms  [Figs.  7(a)  and 
7(b)]  are  identical  to  the  images  reconstructed  with 
the  full-  and  minimal-scan  FBPP  algorithms  [Figs. 
7(c)  and  7(d)]. 

Using  one  of  the  simulated  noisy  complete  and 
minimal-complete  data  sets,  we  again  reconstructed 
the  phantom,  using  the  full-  and  minimal-scan  E-C 
and  FBPP  reconstruction  algorithms.  Figures  8(a) 
and  8(b)  show  the  images  obtained  by  use  of  the 
full-scan  and  the  minimal-scan  E-C  reconstruction 
algorithms,  respectively.  The  images  no  longer  ap¬ 
pear  identical,  and  the  image  reconstructed  with  the 
full-scan  E-C  algorithm  appears  less  noisy  than  the 
image  reconstructed  with  the  minimal-scan  E-C  al¬ 
gorithm.  Figures  8(c)  and  8(d)  show  the  images  ob¬ 
tained  by  use  of  the  f\ill-scan  and  the  minimal-scan 
FBPP  reconstruction  algorithms,  respectively.  Simi¬ 
larly,  the  image  reconstructed  with  the  fiill-scan 
FBPP  algorithm  appears  less  noisy  than  the  image 
reconstructed  with  the  minimal-scan  FBPP  algo¬ 
rithm. 

The  observation  that  the  full-scan  algorithms  gen¬ 
erate  cleaner-looking  images  than  do  the  minimal- 
scan  algorithms  is  not  surprising  and  can  be 
qualitatively  understood  by  examination  of  ways  in 
which  the  redundant  information  inherent  in  the  DT 
data  function  is  utilized.  The  full-scan  FBPP  and 
E-C  algorithms  [with  o>A(vm)  =£  0,  1]  effectively  use 
the  redundant  information  to  reconstruct  two  sepa- 


(c)  (d) 

Fig.  8.  Images  reconstructed  from  noisy  data  by  use  of  (a)  full- 
scan  E-C,  (b)  minimal-scan  E-C,  (c)  full-scan  FBPP,  and  (d) 
minimal-scan  FBPP  reconstruction  algorithms.  The  noisy  data 
were  generated  with  crr  =  a*  =  0.05  in  Eq.  (36). 

rate  images  that  are  averaged  to  form  the  final  image. 
It  has  been  quantitatively  demonstrated  that  this 
effective  averaging  operation  can  result  in  an  unbi¬ 
ased  reduction  of  the  reconstructed  image  vari¬ 
ance.12’14  The  minimal-scan  algorithms,  however, 
utilize  part  of  the  redundant  information  that  is  in¬ 
herent  in  the  data  function  to  reduce  the  angular 
range  over  which  measurements  are  required  for  the 
reconstruction.  The  redundant  information  not 
used  for  this  purpose  can  be  used  to  reduce  the  image 
variance  of  the  reconstructed  image.  Specifically, 
the  complementary  information  contained  in  sub¬ 
spaces  and  <8  of  Fig.  3  is  weighted,  as  described  by 
Eq.  (29),  and  is  subsequently  combined  during  the 
reconstruction  procedure.  However,  unlike  the  full- 
scan  algorithms,  the  minimal-scan  algorithms  cannot 
further  reduce  the  reconstructed  image  variance  by 
exploiting  the  fact  that  subspaces  28  and  9)  contain 
redundant  information. 

Although  the  full-  and  minimal-scan  E-C  algo¬ 
rithms  are  mathematically  equivalent  to  the  full-  and 
minimal-scan  FBPP  algorithms,  we  observed  from 
Fig.  8  that  the  E-C  and  FBPP  algorithms  respond 
differently  to  noise  that  is  present  in  a  discrete  data 
set.  To  confirm  this  observation  quantitatively,  we 
calculated  the  local  image  variances  of  images  recon¬ 
structed,  using  the  different  methods.  Figure  9(a)  is 
a  plot  of  the  local  variance  obtained  from  the 
minimal-scan  FBPP  reconstructed  images  divided  by 
the  local  variance  obtained  from  the  minimal-scan 
E-C  reconstructed  images.  Clearly,  the  ratio  of  the 
variances  is  everywhere  greater  than  1,  quantita¬ 
tively  demonstrating  that  the  minimal-scan  E-C  re¬ 
construction  algorithms  are  less  susceptible  to  the 
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(b) 

Fig,  9,  (a)  Plot  of  the  local  variance  obtained  from  the  full-scan 

FBPP  reconstructed  images  divided  by  the  local  variance  obtained 
from  full-scan  E-C  reconstructed  images,  (b)  Plot  of  the  local 
variance  obtained  from  the  minimal-scan  FBPP  reconstructed  im¬ 
ages  divided  by  the  local  variance  obtained  from  minimal-scan  E-C 
reconstructed  images.  In  both  the  full-  and  the  minimal-scan 
cases,  the  E-C  algorithms  did  a  better  job  of  suppressing  data  noise 
than  did  the  FBPP  algorithms. 


effects  of  data  noise  than  are  the  minimal-scan  FBPP 
reconstruction  algorithms.  Figure  9(b)  is  a  plot  of 
the  local  variance  obtained  from  the  full-scan  FBPP 
reconstructed  images  divided  by  the  local  variance 
obtained  from  the  full-scan  E-C  reconstructed  im¬ 
ages.  The  ratio  of  the  variances  is  everywhere 
greater  than  1,  quantitatively  demonstrating  that 
the  full-scan  E-C  reconstruction  algorithms  are  less 
susceptible  to  the  effects  of  data  noise  than  are  the 
full-scan  FBPP  reconstruction  algorithms.  These 
results  are  consistent  with  those  obtained  previously 
for  full-  and  minimal-scan  DT  utilizing  plane-wave 
illumination.  16>21 

6.  Summary 

In  this  study,  we  revealed  and  examined  the  redun¬ 
dant  information  that  is  inherent  in  the  fan-beam  DT 


data  function.  Such  information  can  be  exploited  to 
reduce  the  reconstructed  image  variance  or  alterna¬ 
tively  to  reduce  the  angular  scanning  requirements  of 
the  fan-beam  DT  experiment.  We  developed  novel 
full-scan  and  minimal-scan  E-C  and  FBPP  recon¬ 
struction  algorithms  for  fan-beam  DT.  The  family  of 
fan-beam  full-scan  E-C  algorithms  operates  by  trans¬ 
forming  (in  2D  Fourier  space)  the  fan-beam  DT  prob¬ 
lem  into  a  2D  parallel-beam  x-ray  CT  problem,  which 
can  be  efficiently  and  stably  inverted  by  use  of  the 
FBPJ  algorithm.  The  family  of  fan-beam  full-scan 
FBPP  algorithms  operates  directly  on  the  modified 
data  function  to  reconstruct  the  image  and  contains 
the  fan-beam  FBPP  algorithm  suggested  by  Dev- 
aney8  as  a  special  member.  Different  members  of 
the  families  of  full-scan  E-C  and  FBPP  algorithms 
are  specified  by  different  choices  of  the  combination 
coefficient  wfe(vm),  which  controls  ways  in  which  the 
redundant  information  in  the  data  function  is  com¬ 
bined.  Reconstruction  algorithms  that  correspond 
to  different  choices  of  c ok(vm)  will  in  general  respond 
differently  to  the  effect  of  noise  and  discrete  sam¬ 
pling.12 

The  fan-beam  minimal-scan  E-C  and  FBPP  algo¬ 
rithms  were  developed  from  the  concept  of  the 
minimal-complete  data  set.  The  minimal-complete 
data  set,  which  is  acquired  by  use  of  view  angles  only 
in  [0,  4>mm]  where  tt  ^  4>min  ^  3tt/2,  contains  all  the 
information  necessary  for  exactly  reconstructing  the 
scattering  object  function.  The  fan-beam  minimal- 
scan  E-C  and  FBPP  algorithms  utilize  a  weighting 
function  w(vm,  c|>)  to  normalize  appropriately  the  par¬ 
tially  redundant  information  inherent  in  the 
minimal-complete  data  set.  Accordingly,  one  can 
form  different  fan-beam  minimal-scan  E-C  and  FBPP 
algorithms  by  specifying  different  choices  for  this 
weighting  function.  Reconstruction  algorithms  that 
correspond  to  different  choices  of  w(vm,  cf>)  will  in 
general  respond  differently  to  the  effect  of  noise  and 
discrete  sampling.15  It  can  be  readily  verified  that, 
under  the  conditions  of  continuous  samplirig  and  in 
the  absence  of  noise,  the  minimal-scan  E-C  and  FBPP 
algorithms  are  exact  and  mathematically  equivalent 
to  their  full-scan  counterparts  that  utilize  measure¬ 
ments  over  the  angular  range  0  <  <j>  <  2tt. 

An  implementation  of  the  fan-beam  full-scan  and 
minimal-scan  algorithms  has  been  presented,  along 
with  numerical  results  obtained  with  noiseless  and 
with  noisy  simulated  data.  It  was  observed  that  the 
full-scan  algorithms  did  a  better  job  of  suppressing 
data  noise  than  did  their  minimal-scan  counterparts. 
We  quantitatively  demonstrated  that  the  full-  and 
minimal-scan  E-C  algorithms  are  less  susceptible  to 
data  noise  and  to  finite  sampling  effects  than  are  the 
full-  and  minimal-scan  FBPP  algorithms,  respec¬ 
tively.  This  result  is  consistent  with  the  observation 
that  the  FBPP-based  algorithms  involve  more- 
complicated  numerical  operations  than  do  the  E-C- 
based  algorithms,  which  may  amplify  the 
propagation  of  noise  and  errors  into  the  reconstructed 
image. 

We  have  assumed  a  2D  imaging  model  in  this 
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study.  Therefore  the  developed  reconstruction  algo¬ 
rithms  may  be  useful  for  applications  in  which  out- 


It  can  readily  be  verified  that  the  four  roots  of  Eq. 
(A2)  are  given  by 


I  Va\l/2[  X4 

vmi=-vm2  =  v^l-^-jj  \1/2{1  -  (1  -  x2)(V/2v02)  +  [1-  X2(l  -  X2)(vo2/v02)]1/2}; 

/.  _ ' 

vm3  Vm4  v0^l  wj  ^y2{1  _  (1  _  x*y(Va*/2vo2)  -  [1  -  X2d  - 


(A3) 

(A4) 


of-plane  scattering  is  not  significant.  The  full-scan 
E-C  and  FBPP  reconstruction  algorithms  can  be  gen¬ 
eralized  readily  to  address  the  three-dimensional  DT 
problem  by  use  of  spherical-wave  sources  and  planar 
measurement  surfaces.  It  remains  unclear  whether 
numerically  stable  versions  of  the  minimal-scan  E-C 
and  FBPP  reconstruction  algorithms  can  be  devel¬ 
oped  for  three-dimensional  imaging  geometries. 

Here  we  have  developed  linear  reconstruction  al¬ 
gorithms  for  fan-beam  DT.  It  was  not  our  intent  to 
address  the  limitations  of  the  Bom  or  Rytov10-11 
weak-scattering  approximation.  The  developed  full- 
and  minimal-scan  algorithms  will,  however,  provide 
a  natural  framework  for  the  incorporation  of  higher- 
order  scattering  perturbation  approximations22-24 
into  the  algorithms.  It  remains  to  be  determined 
whether  minimal-scan  reconstruction  algorithms  can 
be  developed  without  use  of  the  paraxial  approxima¬ 
tion.25  We  intend  to  report  on  the  theoretical  devel¬ 
opment  and  numerical  analysis  of  these  problems  in 
a  forthcoming  publication. 


Appendix  A:  Acronyms  Used 


CT 

DT 

E-C  reconstruction 
algorithm 

FBPJ  reconstruction 
algorithm 

FBPP  reconstruction 
algorithm 

FDP  theorem 


Computed  tomography 
Diffraction  tomography 
Estimate-combination  recon¬ 
struction  algorithm 
Filtered  backprojection  re¬ 
construction  algorithm 
Filtered  backpropagation  re¬ 
construction  algorithm 
Fourier  diffraction  projection 
theorem 


Appendix  B:  Relationships  between  vm  and  va 

From  Eq.  (14)  we  know  that 


For  a  given  va  that  satisfies  0  <  va  <  v0  Vl+l/x^,  we 
would  like  to  find  the  values  of  vm  that  satisfy  Eq. 
(Al).  This  is  equivalent  to  solving  for  the  roots  of  the 
fourth-order  equation: 

a  -  xv(^)‘ + iW  -  2<i  - 

+  v„!(v„!  -  4V,*)  -  0.  (A2) 


Because  when  0  <  va  <  v0Vl  +  1/x2 


21 1/2 


1-X2(1-X2)^2 
v0 


(A5) 


we  observe  that  the  roots  vm3  and  vm4  are  complex 
valued  and  therefore  are  not  physically  meaningful. 
As  expected,  when  x  ->  1,  Eq.  (A3)  reduces  to  the 
known  result12  for  plane-wave  DT  given  by 


/  v02  \1/2 

vmi  =  -vm2  =  vJ  1  -  —~g  I  .  (A6) 
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On  a  Limited- View  Reconstruction  Problem  in  Diffraction 
Tomography 

Xiaochuan  Pan*  and  Mark  A.  Anastasio 


Abstract — Diffraction  tomography  (DT)  is  an  inversion  technique  that 
reconstructs  the  refractive  index  distribution  of  a  scattering  object.  We  pre¬ 
viously  demonstrated  that  by  exploiting  the  redundant  information  in  the 
DT  data,  the  scattering  object  could  be  exactly  reconstructed  using  mea¬ 
surements  taken  over  the  angular  range  [0,  <f>m jn],  where  7r  <  <f>m jn  < 
37t /2.  In  this  paper,  we  reveal  a  relationship  between  the  maximum  scan¬ 
ning  angle  and  image  resolution  when  a  filtered  backpropagation  (FBPP) 
reconstruction  algorithm  is  employed  for  image  reconstruction.  Based  on 
this  observation,  we  develop  short-scan  FBPP  algorithms  that  reconstruct 
a  low-pass  filtered  scattering  object  from  measurements  acquired  over  the 
angular  range  [0, 4>c],  where  Oc  <  <£min. 

Index  Terms — Diffraction  tomography,  limited-view  tomography,  wave- 
field  inversion  techniques. 


I.  Introduction 

In  diffraction  tomography  (DT),  a  semi-transparent  scattering  ob¬ 
ject  is  interrogated  using  a  diffracting  optical  or  acoustical  wavefield 
and  the  scattered  wavefield  around  the  object  is  measured  and  used  to 
reconstruct  the  (low-pass  filtered)  refractive  index  distribution  of  the 
scattering  object.  The  principles  of  DT  have  been  extensively  utilized 
for  developing  optical  and  acoustic  tomographic  imaging  systems.  Re¬ 
cently,  interest  in  DT  within  the  optical  imaging  community  has  in¬ 
creased  because  of  its  potential  application  to  the  diffuse-photon  den¬ 
sity  wave  tomography  [l]-[3]. 

It  was  shown  previously  [4],  [5]  that,  in  two-dimensional  (2-D)  DT 
employing  plane-wave  or  cylindrical- wave  sources  and  the  classical 
scanning  geometry,  one  can  reconstruct  the  scattering  object  from  a 
minimal-scan  data  set  comprised  of  measurements  acquired  over  the 
angular  range  [0,  where  7r  <  <j>m *m  <  37t/2  is  specified  by 

the  measurement  geometry.  In  this  paper,  we  reveal  a  relationship  be¬ 
tween  image  resolution  and  maximum  scan  angle,  based  upon  which 
short-scan  algorithms  can  be  designed  for  reconstructing  a  low-pass 
filtered  scattering  object  from  measurements  acquired  over  the  angular 
range  [0,  <FC],  where  $c  <  <j>m jn.  When  the  scattering  object  is  suf¬ 
ficiently  bandlimited,  it  can  be  exactly  reconstructed  from  the  lim¬ 
ited-view  measurements  in  [0,  $c).  We  present  numerical  examples 
that  confirm  our  theoretical  assertions. 

II.  Background 

Consider  the  classical  scanning  geometry  of  DT  with  a  cylin¬ 
drical  wave  source,  as  shown  in  Fig.  1.  Let  (x,y)  and  (r,  9) 
denote  the  fixed  Cartesian  and  polar  coordinate  systems  and 
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Fig.  1.  The  fan-beam  scanning  geometry  of  2-D  DT.  The  interrogating 
cylindrical  wave  propagates  along  the  77  axis  and  the  scattered  wave  field  is 
measured  along  the  line  77  =  1.  S  (or  D)  denotes  the  distance  between  the 
source  (or  the  detector)  and  the  center  of  rotation. 


(£,  77)  the  rotated  coordinate  system.  These  systems  are  related  by 
x  =  rcosO,  y  =  rsin#,  £  =  xcos <j>  +  ysin<f>  —  rcos(0  —  0),  and 
rj  =  —xsin(l>  +  ycos(j>  =  -rsin(<£  -  9 ).  The  scattering  object,  which 
is  embedded  in  a  lossless  and  homogeneous  background  medium, 
is  illuminated  by  a  monochromatic  cylindrical-wave  ui(£,<j>)  with 
complex  amplitude  Uo  and  wavenumber  k  =  27t^o,  generated  by 
a  line  source  located  at  the  position  77  —  —5  on  the  77  axis.  From 
measurements  of  the  scattered  wavefield  on  the  £  axis  at  different 
view  angles  </>,  one  seeks  to  reconstruct  the  scattering  object  function 
a(f),  which  is  related  to  the  refractive  index  distribution  n(r)  within 
the  scattering  object  by  a(r)  =  n2(r)  -  1. 

Let  u(£,  <f> )  and  us  (£,  <j>)  —  u(£,  <j>)  -  Ui(£,  <j>)  denote  the  total  and 
scattered  wavefields  measured  along  the  line  rj  =  D  oriented  at  angle 
<t>,  as  shown,  in  Fig.  1.  For  the  sake  of  convenience,  we  introduce  a 
modified  data  function  M(i/m,  <j>)  that  can  be  obtained  readily  from 
the  scattered  wavefield  and  is  defined  as 

M(vm,  <j>)  =  -^2  "'exp  [— j2ir(j/  -  u0)D]  (1) 

where  *  =  S/JsTD ),  v'  =  0/g  -  ^/x2,  and  ^m{h(f)}.  = 
1  /(2tt)  f  ^  h(£)e  d£.  The  special  case  of  plane-wave  illumi¬ 

nation  (S  — ►  00)  corresponds  to  x  =  1-  Under  the  Bom  and  paraxial 
approximations,  Devaney  derived  the  fan-beam  Fourier  diffraction  pro¬ 
jection  (FDP)  theorem  [6],  which  relates  a(r)  to  the  modified  data 
function  by 


/oo  r  00 

/  a(f> 

-OO  7  —  00 

'  exp  j-j27T  ^~£  -  (1/  -  v0)  7}  |  dr, 


if  l^m |  <X^o 

=0  ifkm|>X^o- 


(2) 
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The  FDP  theorem  can  also  be  derived  by  employing  the  Rytov 
approximation.  In  this  case,  (2)  remains  unchanged  and  only  (1)  needs 
to  be  appropriately  redefined  [6].  When  <f>  is  varied  from  0  to  27t,  the 
FDP  theorem  specifies  a  circular  disk  of  (double)  coverage  centered 
at  the  origin  with  radius  vo  y/l  T  1/x2  *n  the  2-D  Fourier  space  of 
a(r).  Conventional  full-scan  reconstruction  algorithms,  such  as  the 
well  known  filtered  backpropagation  (FBPP)  algorithm  [7],  utilize  this 
Fourier  space  coverage  for  reconstructing  a(r).  We  will  refer  to  such 
a  (low-pass  filtered)  reconstructed  a(r)  as  the  “exact”  image. 

According  to  the  fan-beam  FDP  theorem  in  (2),  the  modified  data 
function  M(vm,4>)  satisfies  the  consistency  condition  [4] 


M  (Vm,  4>)  =  M  {-VmA  +  7T  -  2a)  (3) 

where  sin  a  =  sgn(t'm)  [((i/  -  ^o)2/(^m/x2)  4-  (?'  -  ^o)2)] 1/2- 
Using  (3),  one  can  show  [4],  [5]  that  the  minimal-scan  data  acquired 
in  the  angular  range  [0,  ^m;n]  specifies  a  circular  disk  (with  radius 
v o  y/l  4-  1/x2)  of  coverage  in  the  Fourier  space  of  a(f),  where 


(t> min  =  7T  4*  2 8  and 


(4) 


III.  A  Limited-View  Reconstruction  Problem  for  2-D  DT 


-XVQ  -pc  0  VC  XVO 

Vm 

Fig.  2.  The  full-scan  data  space  Wc  =  A\J  BDC\JV  contains  data  in  the 
angular  range  [0,  27r].  The  subspace  A  U  B  U  C  in  [0,  4>c]  contains  all  of  the 
information  ncessary  for  exact  reconstruction  of  the  scattering  object  function 
whose  2-D  Fourier  transform  is  bandlimited  to  a  disk  of  radius  Rc(ue). 


We  focus  now  on  a  limited- view  problem,  in  which  data  are  acquired 
only  over  the  angular  range  [0,  $c],  where  7r  <  $c  <  </>min.  In  this 
situation,  it  is  well  known  that  the  exact  image  cannot,  in  general,  be  re¬ 
constructed  [8].  However,  we  demonstrate  that  algorithms  can  be  devel¬ 
oped  for  reconstructing  a  low-pass  filtered  approximation  of  the  exact 
image.  Consider  a  scattering  object  ac(r)  whose  2-D  Fourier  transform 
Ac{v)  is  bandlimited  to  a  disk  of  radius  Rc  centered  at  the  origin,  where 


and  0  <  vc  <  xvo  ■  Then,  according  to  the  fan-beam  FDP  theorem 
in  (4),  the  modified  data  function  M(vm,  <f>)  is  nonzero  only  for 
\vm |  <  vc.  The  data  space  Wc  =  [|*/m|  <  vci0  <  <j>  <  2i r],  in 
which  the  modified  data  function  is  defined,  can  be 

divided  into  the  four  subspaces  A,  By  C,  and  X>,  as  shown  in 
Fig.  2,  where  A  ~  [|^m|  <  vc, ,  0  <  <j>  <  2a(vm)  +  2a(i/c)J, 

B  =  [|^m|  <  vc ,  2a(i/m)  +  2a(z/c)  <  <j>  <  tt  +  2a(i/m)j, 

C  =  [\vm\  <vc,TT  +  2 a(vm)  <4><  ^c],  and 

V  =  [\vm\  <1 'c,$c  <<j>  <  2tt].  The  value  of  $c  is  determined  by 

$c  =  7r  +  2a  (vc)  .  (6) 

Using  (3),  it  can  be  verified  that  information  of  M(vm,  <t>)  in  subspace 
A  is  redundant  to  that  of  in  subspace  C.  Similarly,  infor¬ 
mation  of  (j>)  in  subspace  B  is  redundant  to  that  of  4>) 

in  subspace  V.  Therefore,  in  principle,  the  modified  data  function 
AT (^m,  d>)  is  completely  specified  by  its  values  in  the  subspaces  A 
and  B.  However,  because  the  boundary  between  the  subspaces  B  and 
C  is  a  nonlinear  function  of  Vm  and  <j>  and  because  each  horizontal 
line  in  Wc  corresponds  to  a  measurement  acquired  at  a  particular 
angle  <j>,  the  information  in  subspaces  B  and  C  cannot  in  practice  be 
determined  independently  of  each  other.  Consequently,  in  order  to 
determine  M(z/m,^)  in  subspaces  A  and  B,  it  is  necessary  to  scan 
the  union  AU  BUC  =  [|i/m|  <  vc,  0  <  <j>  <  $c].  This  observation 
can  also  be  understood  by  examining  the  2-D  Fourier  space  coverage 
of  a(r)  that  is  obtained  by  varying  the  scanning  angle  from  0  to  $c. 
As  shown  in  Fig.  3,  although  the  disk  of  Fourier  space  coverage  with 


Fig.  3 .  The  2-D  Fourier  space  coverage  of  the  scattering  object  that  is  obtained 
by  varying  <j>  from  0  to  4>c.  The  disk  of  Fourier  space  coverage  with  radius 
R  —  v o  yj\  4  1/x2  is  incomplete,  with  the  shaded  region  denoting  the  missing 
data.  However,  the  coverage  corresponding  to  the  disk  of  radius  Rc ,  which  is 
defined  by  (5),  is  completely  specified.  In  generating  the  figure,  $c  =  230° 
and  x  =  1  were  utilized. 

radius  R  =  vo  y/l  +  1/x2  is  incomplete,  the  coverage  corresponding 
to  the  disk  of  radius  Rc{vc)  is  completely  specified.  Therefore,  in 
order  to  exactly  reconstruct  ac(r)  whose  2-D  Fourier  transform 
Ac(v)  is  bandlimited  to  a  disk  of  radius  Rc(vc),  only  measurements 
corresponding  to  view  angles  in  [0,  4>c]  are  required.  Alternatively, 
for  an  arbitrary  scattering  object  a(r)  and  specified  $c  >  7r,  one 
can  readily  reconstruct  ac(f),  which  is  a  low-pass  filtered  version  of 
a(r)  whose  2-D  Fourier  transform  is  bandlimited  to  the  disk  of  radius 
Rc{vc)>  where  the  value  of  the  data  cutoff  frequency  0  <  vc  <  x^o 
is  determined  by  (6). 

A  plot  of  $c  versus  vc/v 0  for  plane-  and  cylindrical- wave  illumina¬ 
tion  is  shown  in  Fig.  4.  As  expected,  the  maximum  scanning  angle  <£c 
is  a  monotonically  increasing  function  of  the  data  cutoff  frequency  vc. 
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The  nonlinear  shape  of  the  curves  indicates  the  scanning  angle  can  be 
reduced  from  (e.g.,  by  30°)  with  little  loss  of  resolution  in  the  re¬ 
constructed  image.  Also,  the  fact  that  the  plane-wave  (x  =  1)  curve  is 
everywhere  higher  than  the  cylindrical-wave  (x  =  0.5)  curve  reflects 
the  fact  that  the  angular  scanning  requirements  of  plane-beam  DT  are 
more  restrictive  than  for  fan-beam  DT  [5]. 


IV.  Short-Scan  Reconstruction  Algorithms  for  Limited-View 

DT 

Although  the  data  AU  B  U  C  in  Fig.  2  contains  all  of  the  informa¬ 
tion  necessary  for  reconstruction  of  ac(f),  subspaces  A  and  C  contain 
redundant  information  that  needs  to  be  properly  normalized  in  the  re¬ 
construction  process.  This  can  be  achieved  by  introducing  a  weighted 
modified  data  function  as  [4],  [5] 

M*  {Vm,  <t>)  =  l v(Vm,<t>)M  (I'm,  <t>)  (7) 

where  w(i/m,  <j>)  satisfies 

w  (i/m,  4>)  +  w  (  I'm ,  4>  +  7r  —  2a)  =  1  (8a) 

everywhere  in  the  data  space  Wc 

w  (I'm,  6)  (8b) 

in  subspace  B  and 

w{Vm,<)>)  =  0  (8C) 


in  the  subspace  {D  U  [\vm  |  >  ^,0  <  <j>  <  2tt]}.  The  image  ac(r)  can 
be  reconstructed  using  a  short-scan  FBPP  (SS-FBPP)  reconstruction 
algorithm  given  by 


a(m\r,9)=  f*  F  ^\vm\ 

J  <jx=0  J I/m=5 —  Vc  V 


X  exp  \j2ir  sgn (vm)u  ^  +(j'/-  i/q)  rcos(<£-  a -  9) 


which  reduces  to  the  full-scan  fan-beam  FBPP  algorithm  [5]  when 
$c  =  27r  and  w(um,<t>)  =s  1/2.  Note  that  different  choice  for 
w(vm,(t>)  that  satisfy  (8),  in  effect,  specify  different  SS-FBPP 
algorithms. 


V.  Numerical  Results 

To  validate  the  theoretical  results  above,  we  considered  a  numerical 
phantom  containing  two  elliptical  disks  whose  2-D  Fourier  transform 
was  approximately  bandlimited  to  a  disk  of  radius  Rc{vc  —  0.451/0 ) 
[see  (5)].  Data  sets  of  simulated  scattered  fields  were  generated  using 
the  plane- wave  FDP  theorem  (i.e.,  x  =  1)  and  using  various  values  for 
$c.  We  reconstructed  images,  which  are  shown  in  Fig.  5,  from  these 
data  sets  using  the  conventional  FBPP  and  SS-FBPP  algorithms.  The 
SS-FBPP  algorithm  was  specified  by  a  weighting  function  w(vm,  <j>) 
that  took  on  the  values  1/2, 1 , 1/2,  and  0,  in  the  data  subspaces  A ,  B ,  C, 
and  V ,  respectively.  Fig.  5(a)  shows  images  reconstructed  by  use  of 
the  FBPP  algorithm  (left)  and  SS-FBPP  algorithm  (right),  using  data 
sets  corresponding  to  $c  =  2n  and  <f>c  =  <£min(x  =  1)  =  37r/2, 
respectively.  It  is  observed  that  both  images  appear  virtually  identical, 
reflecting  the  fact  that  both  of  these  data  sets  contain  the  complete  in¬ 
formation  about  the  scattering  object. 

Fig.  5(b)  shows  images  reconstructed  by  use  of  the  FBPP  algorithm 
(left)  and  SS-FBPP  algorithm  (right),  using  a  data  set  corresponding 
to  $c  =  7r  4-  2a(0.45i/oJ(«  207°).  Clearly,  the  image  reconstructed 

using  the  FBPP  algorithm  is  distorted  and  contains  artifacts.  However, 
the  image  reconstructed  by  use  of  the  SS-FBPP  algorithm  appears  cor¬ 
rect  and  virtually  identical  to  the  images  shown  in  Fig.  5(a).  This  con¬ 
firms  our  assertion  that  the  SS-FBPP  algorithms,  which  utilize  in  the 
angular  range  [0,  $c],  can  exactly  reconstruct  a  scattering  object  ac(r) 
whose  2-D  Fourier  transform  Ac(i /)  is  bandlimited  to  a  disk  of  radius 
flc(i'c),  where  vc  and  $c  are  related  by  (6). 

Fig.  5(c)  shows  images  reconstructed  by  use  of  the  FBPP  algorithm 
(left)  and  SS-FBPP  algorithm  (right),  using  a  data  set  corresponding  to 
<$>c  =  x  -f  2a;(0.25i'o)(«  195°).  Note  that  because  the  2-D  Fourier 
transform  of  a(r)  has  support  on  the  disk  or  radius  Rc( vc  =  0.45i/o)> 
the  measurements  in  the  angular  range  [0,  3>c  =  195°]  do  not  com¬ 
pletely  specify  the  scattering  object  (i.e.,  the  disk  of  coverage  in  2-D 
Fourier  space  with  radius  Rc(yc  =  CUS^o)  will  not  be  completely 
filled  in.)  As  expected,  the  image  reconstructed  using  the  FBPP  al¬ 
gorithm  is  blurred,  distorted  and  contains  artifacts.  The  image  recon¬ 
structed  using  the  SS-FBPP  algorithm  also  appears  blurred,  but  does 
not  contain  any  noticeable  distortions  or  artifacts.  This  confirms  our 
assertion  that,  when  the  2-D  Fourier  transform  of  a  scattering  object 
a(f)  is  not  bandlimited  to  a  disk  of  radius  Rc(vc),  the  SS-FBPP  algo¬ 
rithms  that  utilize  the  measurements  corresponding  to  view  angles  in 
[0,  $c]  (where  v0  and  <$>c  are  related  by  (6))  can  reconstruct  a  low-pass 
filtered  version  of  a(r)  whose  2-D  Fourier  transform  is  bandlimited  to 
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(C) 


Fig.  5.  Images  reconstructed  using  the  FBPP  and  SS-FBPP  algorithms  for 
various  simulated  data  sets.  See  the  text  for  a  detailed  description. 


the  disk  of  radius  Rc{vc ).  In  this  particular  example,  the  2-D  Fourier 
transform  of  the  image  reconstructed  by  use  of  the  SS-FBPP  algorithm 
is  bandlimited  to  a  disk  of  radius  i?c(.25/.45  vc ). 

VI.  Summary 

We  demonstrated  previously  [4],  [5]  that  in  2-D  DT  employing 
plane-wave  or  cylindrical-wave  sources,  one  can  exactly  reconstruct 
the  scattering  object  from  a  minimal-scan  data  set  acquired  using  view 
angles  only  in  [0,  0min],  where  it  <  (j>m] in  <  3ir/2  is  a  specified 
function  of  the  measurement  geometry.  In  this  study,  we  have  demon¬ 
strated  that  when  measurements  are  available  only  for  view  angles  in 
[0,  4>c],  where  tt  <  $c  <  a  simple  relationship  exists  between 
the  maximum  scanning  angle  4>c  and  the  image  resolution  when  a 
FBPP  algorithm  is  employed  to  reconstruct  the  image.  By  properly 
weighting  the  measurement  data,  a  low-pass  filtered  approximation 
of  the  scattering  object  that  is  free  of  conspicuous  artifacts  can  be 
obtained  from  the  measurements  corresponding  to  view  angles  in 
[0,  $c].  When  the  scattering  object  is  sufficiently  bandlimited,  it 
can  be  exactly  reconstructed.  This  observation  is  practically  useful, 
because  it  provides  a  convenient  mechanism  for  regularizing  the 
severely  ill-posed  limited- view  DT  reconstruction  problem;  when  the 
maximum  scanning  angle  4>c  is  greater  than  7r,  a  stable  reconstruction 
can  always  be  performed  by  sacrificing  spatial  resolution  in  the  recon¬ 
structed  image.  It  can  be  demonstrated  that  the  statistical  properties 
of  the  SS-FBPP  algorithms  are  qualitatively  similar  to  those  of  the 
minimal-scan  FBPP  reconstruction  algorithms  investigated  previously 
[5].  In  the  limited- view  radon  transform  inversion  problem  [9],  an 
analogous  regularization  mechanism  does  not  exist  and  some  sort  of  a 
priori  information  regarding  the  object  function  is  generally  required 
to  effectively  regularize  the  problem. 

Because  we  have  assumed  a  2-D  imaging  geometry  in  this  study,  the 
developed  SS-FBPP  reconstruction  algorithms  may  be  useful  for  appli¬ 
cations  in  which  out-of-plane  scattering  is  not  significant.  In  diffuse- 
photon  density  wave  tomography,  the  wavenumber  is  complex-valued 
and  the  FDP  theorem  describes  a  mapping  between  the  data  function 
and  a  set  of  complex -valued  frequencies  of  the  scattering  object  func¬ 
tion’s  Fourier  transform.  The  extension  of  the  concepts  and  techniques 
introduced  in  this  correspondence  to  the  case  where  the  wavenumber 
is  complex-valued  and  to  the  three-dimensional  reconstruction  problem 
represent  important  topics  for  future  research. 
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ABSTRACT:  It  is  widely  believed  that  measurements  from  a  full 
angular  range  of  are  generally  required  to  exactly  reconstruct  a 
complex-valued  refractive  index  distribution  in  diffraction  tomography 
(DT).  In  this  work,  we  developed  a  new  class  of  minimal-scan  recon¬ 
struction  algorithms  for  DT  that  utilizes  measurements  only  over  the 
angular  range  0  <  <f>  <  3W2  to  perform  an  exact  reconstruction. 
These  algorithms,  referred  to  as  minimal-scan  estimate-combination 
(MS-E-C)  reconstruction  algorithms,  effectively  operate  by  transform¬ 
ing  the  DT  reconstruction  problem  into  a  conventional  x-ray  CT  re¬ 
construction  problem  that  requires  inversion  of  the  Radon  transform. 
We  performed  computer  simulations  to  compare  the  noise  and  nu¬ 
merical  properties  of  the  MS-E-C  algorithms  against  existing  filtered 
backpropagation-based  algorithms.  ©  2002  Wiley  Periodicals,  Inc.  Int  J 
Imaging  Syst  Technol,  12,  84-91,  2002;  Published  online  in  Wiley  InterScience 
(www.interscience.wiley.com}.  DOI  10.1002/ima.10014 

Key  words:  topographic  reconstruction;  diffraction  tomography; 
wavefield  inversion 


i.  INTRODUCTION 

In  diffraction  tomography  (DT),  a  scattering  object  is  interrogated 
using  a  diffracting  acoustical  or  electromagnetic  wavefield,  and  the 
scattered  wavefield  around  the  object  is  measured  and  used  to 
reconstruct  the  refractive  index  distribution  of  the  scattering  object. 
There  are  numerous  potential  applications  of  DT  that  can  be  found 
in  various  scientific  fields  (Andre  et  al.,  1995;  Tabbara  et  al.,  1988; 
Mueller  et  ah,  1979;  Kino,  1979;  Devaney,  1984;  Robinson,  1984). 
Recently,  there  has  also  been  considerable  interest  in  using  DT  to 
perform  coherent  x-ray  imaging  using  third-generation  synchrotron 
sources  (Cheng  and  Han,  2001).  Unlike  the  x-rays  used  in  computed 
tomography  (CT)  that  travel  along  straight  lines,  the  radiation  em¬ 
ployed  in  DT  has  to  be  treated  in  terms  of  wavefronts  and  fields 
scattered  by  inhomogeneities  in  the  object.  In  DT,  the  interaction 
between  the  incident  wavefield  and  the  object  medium  is  governed 
by  the  inhomogeneous  Helmholtz  equation.  Using  a  weak-scattering 
approximation,  the  inhomogeneous  equation  can  be  analytically 
solved  (Wolf,  1969;  Mueller  et  al.,  1979)  to  obtain  a  linear  relation¬ 
ship  between  the  scattered  field  and  the  refractive  index  distribution. 
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This  relationship  has  been  used  to  develop  DT  reconstruction  algo¬ 
rithms  such  as  the  well-known  filtered  backpropagation  (FBPP) 
algorithm  (Devaney,  1982),  which  is  a  generalization  of  the  filtered 
backprojection  (FBPJ)  algorithm  of  x-ray  CT. 

It  is  widely  believed  that  measurements  from  a  full-angular 
range  of  27r  around  the  scattering  object  are  generally  required  to 
exactly  reconstruct  a  complex-valued  refractive  index  distribu¬ 
tion  (Devaney,  1982).  However,  we  have  recently  revealed  that 
one  needs  measurements  only  over  the  angular  range  0  <  <f>  < 
3tt/2  to  perform  an  exact  reconstruction,  and  we  developed 
minimal-scan  filtered  backpropagation  (MS-FBPP)  algorithms  to 
achieve  this  (Pan  and  Anastasio,  1999).  A  useful  characteristic  of 
the  MS-FBPP  algorithms  is  their  ability  to  decrease  the  data 
acquisition  time  by  at  least  25%  over  conventional  (full-scan) 
algorithms.  They  can  also  reduce  artifacts  due  to  movement  in  or 
temporal  fluctuations  of  the  scattering  object.  Furthermore,  in 
certain  practical  situations,  it  may  be  impossible  to  acquire 
measurements  over  a  full  2 tt  angular  range. 

A  new  class  of  reconstruction  algorithms  has  recently  been 
developed  for  full-scan  DT  (in  other  words,  DT  employing  mea¬ 
surements  over  a  full  27t  angular  range).  These  algorithms,  referred 
to  as  estimate-combination  (E-C)  reconstruction  algorithms  (Pan, 
1998;  Anastasio  and  Pan,  2000b;  Anastasio  and  Pan,  2000a),  effec¬ 
tively  operate  by  transforming  the  DT  reconstruction  problem  into  a 
conventional  x-ray  CT  reconstruction  problem  that  can  be  efficiently 
solved  using  the  filtered  backprojection  (FBPJ)  algorithm.  The  E-C 
reconstruction  algorithms  are  more  computationally  efficient  than 
the  FBPP  algorithm,  and  also  provide  a  flexible  framework  for 
imposing  unbiased  regularization. 

Because  the  E-C  reconstruction  algorithms  involve  a  Fourier 
series  expansion  of  the  data  function  that  is  acquired  over  the 
angular  range  0  <  <  2tt,  they  cannot  be  applied  directly  to  the 

minimal-scan  problem  where  measurements  are  only  acquired  over 
the  angular  range  0  <  <f>  <  37t/2.  Because  of  the  potential  advan¬ 
tages  of  the  E-C  reconstruction  algorithms,  it  is  important  to  gen¬ 
eralize  them  to  the  minimal-scan  situation.  In  this  work,  we  devel¬ 
oped  rninimal-scan  E-C  (MS-E-C)  reconstruction  algorithms  for  DT. 
We  performed  computer  simulations  to  compare  the  noise  and 
numerical  properties  of  the  MS-E-C  and  MS-FBPP  algorithms.  Our 
results  quantitatively  demonstrate  that  the  MS-E-C  algorithms  pos- 
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Figure  1 .  The  classical  scanning  geometry  of  2D  DT.  The  insonifying 
plane  wave  propagates  along  the  rj  axis,  and  the  scattered  wave  field 
is  measured  along  the  line  77  =  I.  Full-scan  and  minimal-scan  data  sets 
are  obtained  by  varying  the  measurement  angle  <f>  between  0  and  2 rr 
or  between  0  and  3tt/2,  respectively. 

sess  statistical  and  numerical  properties  superior  to  those  of  the 
MS-FBPP  algorithms. 

II.  BACKGROUND 

A.  The  Fourier  Diffraction  Projection  Theorem.  In  two-di¬ 
mensional  (2D)  DT  employing  the  classical  scanning  configura¬ 
tion,  as  shown  in  Figure  1,  the  scattering  object  is  illuminated  by 
monochromatic  plane-wave  radiation  of  frequency  v0 ,  and  the 
transmitted  wavefield  is  measured  along  the  £  axis  oriented  at  a 
measurement  angle  4>,  at  a  distance  1 7  =  /  from  the  origin.  From 
measurements  of  the  scattered  wavefield  obtained  at  various 
angles  <f>,  one  seeks  to  reconstruct  the  scattering  object  function 
a(r ),  which  is  related  to  the  refractive  index  distribution  n(r,  0)  by 
a(r,  0)  =  n2(r,  0)  -  1. 

At  a  measurement  angle  </>,  the  scattered  data  are  measured  along 
the  line  7]  =  l,  as  shown  in  figure  1.  Let  Us(vm,  4>)  to  denote  the 
ID  Fourier  transform  of  the  measured  scattered  data  with  respect  to 
For  convenience,  we  define  a  modified  ID  Fourier  transform  of 
the  scattered  data  as 

<*>)  =  UAvm  <£>)  2^u0  e~jWI'  (1) 

where  v*  -  Vvq  -  v2  and  |vj  <  v0.  The  quantities  on  the 
right-hand  side  of  equation  1  are  known  or  can  be  measured. 
Therefore,  we  will  treat  M{vnv  (f>)  as  a  measurable  data  function. 
Under  the  Bom  approximation  (Mueller  et  al.,  1979),  the  Fourier 
diffraction  projection  (FDP)  theorem  (Mueller  et  al.,  1979)  can  be 
derived,  which  is  mathematically  stated  as 


0  if  \vm\  >  v0, 

(2) 


where  the  polar  coordinates  (r,  0)  and  the  rotated  coordinates  (£,  17) 
are  related  through  £  =  r  cos(c£  —  0)  and  77  ~  -r  sin (<£  -  0).  The 
FDP  theorem  indicates  that  M(vm,  <f>)  provides  the  values  of  the  2D 
Fourier  transform  of  a(r)  along  the  semi-circular  arc  AOB  of  radius 
v0,  as  shown  in  figure  2. 

B.  Minimal-Scan  Filtered  Backpropagation  Algorithms. 

The  widely  used  filtered  backpropagation  (FBPP)  algorithm  (Dev- 
aney,  1982)  is  mathematically  expressed  as 

|  f  2ir  f  »o  y 

a(r)  =  -  ^  \vm\M(v„„  ft,*™#*™  dv,„  d<f>, 

J  <f)=0  J  vm=-VQ 

(3) 

where  v ^  =  j(Vv \  -  vf„  -  v0).  When  v0  —>  the  FBPP 
algorithm  reduces  to  the  filtered  backprojection  (FBPJ)  algorithm  of 
x-ray  CT.  The  FBPP  algorithm  generally  requires  full  knowledge  of 
M(vm,  <f>)  in  the  data  space  W  =  [\vm\  <  v0,  0  <  <f>  <  27t],  for 
exact  reconstruction  of  the  generally  complex-valued  object  func¬ 
tion.  We  will  refer  to  such  full  knowledge  of  M(v,„,  4>)  as  a  full-scan 
data  set .  The  full-scan  data  space  Tf  can  be  decomposed  into  four 
subspaces,  si,  2fc,  c€,  and  3),  where  si  -  [\vm\  <  v0 ,  ,  0  ^  4>  ^ 
2a  +  7t/2],  2&  =  [|j/m|  ^  v0 ,  tt/2  4-  2a  <  <#>  <  77  +  2a),  % 
=  [\vm\  ^  v0,  tt  +  2a  <  <f>  <  3tt/2]  and  9)  =  [\vj  <  v0 , 
37t/2  <  $  —  2tt],  where  a  ~  sgn(vm)arcsinVv2  -  v2J(2v0).  A 
schematic  of  this  partitioning  of  the  data  space  is  given  in  figure  3. 
Using  the  FDP  theorem,  it  can  be  shown  (Pan  and  Kak,  1983)  that 
M(vm,  4>)  =  M(  —  vm,  4>  +  ir  -  2a).  Therefore,  the  information 
contained  in  subspace  si  is  redundant  to  that  contained  in  subspace 
%,  and  the  information  contained  in  subspace  2&  is  redundant  to  that 
contained  in  subspace  2).  We  have  demonstrated  (Pan  and  Anasta- 
sio,  1999)  that  it  is  possible  to  exactly  reconstruct  the  object  function 
using  only  knowledge  of  M(vm ,  4>)  in  the  subspace  M  —  si  U  2& 


Figure  2.  The  FDP  theorem  states  that  M{vm,  <f>)  is  equal  to  the  2D 
Fourier  transform  of  a(r)  along  a  semi-circle  AOB  that  has  a  radius  of 
v0  and  is  centered  at  =  v0. 
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Figure  3.  The  subspaces  d,  %  <6  and  2  in  the  complete  DT  data 
space. 


U  *€  -  [|v„J  ^  v0,  0  <  </>  <  3tt/2],  which  we  refer  to  as  a  a 
minimal-scan  data  set. 

Although  the  minimal-scan  data  set  it  contains  all  of  the  infor¬ 
mation  necessary  for  exact  reconstruction  of  the  scattering  object 
function,  the  redundant  information  contained  in  the  subspaces  si 
and  %  needs  to  be  properly  normalized  in  the  reconstruction  process 
(Pan  and  Anastasio,  1999).  The  MS-FBPP  algorithms  operate  by 
first  normalizing  such  partially  redundant  information  by  generating 
an  appropriately  weighted  minimal-scan  data  set  M’(vn}y  <j>)  and 
subsequently  using  the  FBPP  algorithm  described  by  equation  3 
(scaled  by  a  factor  of  2),  to  exactly  reconstruct  the  image.  The 
weighted  minimal-scan  data  set  is  given  by  (Pan  and  Anastasio, 
1999) 

M’(vm  =  w(vm%  </>),  (4) 

where  w{vm,  4>)  is  a  function  of  vm  and  </>,  which  satisfies 

w(^„,  4)  +  w{-vm9  t/>  +  rr  -  2a)  =  1  (5a) 


in  subspace  2.  Although  the  forms  of  vr(  vm,  $)  in  subspaces  and 
2)  are  completely  specified  by  equations  5b  and  5c,  respectively,  the 
explicit  forms  of  w(v„r  <£)  in  subspaces  d  and  <6  are  unspecified  for 
the  moment.  In  principle,  one  can  choose  different  w(vnr  cf>)  in 
subspaces  d  and  %  as  long  as  these  w(vnr  <J>)  satisfy  equation  5a. 

111.  MINIMAL-SCAN  ESTIMATE-COMBINATION 
ALGORITHMS 

The  previously  derived  (full-scan)  E-C  reconstruction  algorithms  are 
more  computationally  efficient  than  the  (full-scan)  FBPP  reconstruc¬ 
tion  algorithms,  which  involve  a  depth-dependent  filtering  operation 
(backpropagation).  Accordingly,  we  expect  that  the  MS-FBPP  algo¬ 
rithms,  which  use  the  FBPP  algorithm  to  reconstruct  the  final  image 
from  the  weighted  data  function  M’(vm ,  <£),  will  also  be  less 
computationally  efficient  than  the  E-C  reconstruction  algorithms. 
Because  they  will  involve  fewer  and  less  complicated  numerical 
operations,  we  also  expect  that  the  MS-E-C  algorithms  will  be 
minimize  the  propagation  of  data  noise  and  errors  as  compared  to 
the  MS-FBPP  algorithms.  The  full-scan  E-C  reconstruction  algo¬ 
rithms  (Pan,  1998;  Anastasio  and  Pan,  2000b)  involve  a  Fourier 
series  expansion  of  the  data  function  M{vm ,  </>),  which  requires 
knowledge  of  M(vm,  </>)  over  an  angular  range  of  2 tt,  and  therefore 
can  not  be  directly  applied  to  the  minimal-scan  data  set  containing 
only  measurements  in  the  range  0  <  (f>  <  3tt/3.  Below,  we  develop 
minimal-scan  E-C  (MS-E-C)  reconstruction  algorithms  that  can  be 
directly  applied  to  the  minimal-scan  data  set. 

A.  The  Radon  Transform.  Let  p(f,  <J>)  and  Pk(  vet)  denote  the 
Radon  transform  of  a(r,  0)  and  its  2D  Fourier  transform,  respec¬ 
tively.  (Here,  the  2D  Fourier  transform  is  actually  a  ID  Fourier 
transform  with  respect  to  £  and  a  1 D  Fourier  series  with  respect  to 
<M  Front  knowledge  of  /?(£,  <£),  or  equivalently,  Pk{va ),  one  can 
reconstruct  a{r ,  0)  by  use  of  the  computationally  efficient  and 
numerically  stable  FBPJ  algorithm,  which  is  given  by 

]  f  2ir  f*  * 

a(r,6)  =  -  ^  Pk(va)eik*\v„\ei,,‘K°"*-,)  dv„  d<f>. 

j  0  J  -*  A— * 

(6a) 

For  theoretical  convenience,  the  FBPJ  algorithm  can  also  be  ex¬ 
pressed  as 


tf(r,  0)  -  2 7T  2  /  Pk(v()eikeJk(2irvar)vadvtn  (6b) 

*=-*  J^o 

where  Jk(  *  )  is  a  Bessel  function  of  the  first  kind.  The  MS-E-C 
algorithms  will  operate  by  estimating  Pk(va)y  or  equivalently  p(£, 
<£),  from  the  minimal-scan  data  set,  and  using  the  FBPJ  algorithm  to 
reconstruct  the  final  image  a{ry  0). 


in  complete  data  space  °lf, 


h*(vw.  </))  =  1 


in  subspace  2ft,  and 


«’( K-  4>)  =  0 


(5b) 


(5c) 


B.  Derivation  of  the  MS-E-C  Algorithms.  Consider  a  given 
weighting  function  \v(  vm,  </>)  in  equation  4.  The  corresponding 
MS-FBPP  algorithm  can  be  expressed  as  (Pan  and  Anastasio,  1999) 


a(r,  0)  =  : 


1 


7  k,l <t>)ei2n^‘“ridvmd4>. 


(7) 
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Let  M'k{vm)  denote  the  Fourier  series  expansion  of  M’{vm ,  $).  One 
can  re-express  equation  7  as 

a(r,  0) 


^  M'k{vn)^dvmd<\>. 


(8) 


Using  the  definition  7(1^,)  =  eyof  and  separating  the  contribution  to 
the  integral  from  positive  and  negative  v,„ ,  equation  8  can  be 
re-written  as 


a(r,  6)  =  i  f  f  ^  VmgP-^-^M  £  m;.(0  dv„,  d<t> 

^  6=0  ^  v,„-0 


x  2  M'k(vm)  dv,„  d<t>.  (9) 
*=-* 


Changing  v,„  to  —  in  the  second  term  in  equation  9  and  grouping 
^-dependent  terms  yields 


a(r,  0)  =  j 


ej2irvl„£+2iTvtt7)+jk<t>  I 


X  2  M'k(vn)  dv,„ 

J t=-2C 


X 


g  -;2irvH,f+2ir^Tf+yA-«f> 


2  M[(-vm)dvm. 

J t=-* 


(10) 


and  therefore  the  Radon  transform  of  the  scattering  object  function 
a(r,  0)  can  be  estimated  from  the  appropriately  weighted  minimal- 
scan  data  set.  The  use  of  equation  12  to  estimate  Pk(va)  coupled  with 
the  2D  FBPJ  algorithm  to  reconstruct  a(r ,  6)  is  referred  to  as  a 
MS-E-C  reconstruction  algorithm.  In  practice,  the  FBPJ  algorithm 
described  by  equation  6a  requires  knowledge  of  Pk(va)  for  evenly 
spaced  values  of  va  spanning  the  range  —  V2  v0  ^  va  ^  V2  v0 
in  order  to  be  efficiently  implemented  using  the  fast  Fourier  trans¬ 
form  (FFT).  In  this  case,  the  consistency  condition  (Deans,  1983) 
Pk(va)  =  (-1)*P*(- vfl)  can  be  employed  to  obtain  values  of 
Pk( va)  for  negative  va. 

IV.  NUMERICAL  SIMULATIONS 

We  performed  simulation  studies  to  evaluate  and  compare  the  nu¬ 
merical  and  statistical  properties  of  images  obtained  by  use  of  the 
MS-E-C  and  MS-FBPP  reconstruction  algorithms. 

A.  Measurement  Data.  We  investigated  the  statistical  proper¬ 
ties  of  the  reconstruction  algorithms  under  near-ideal  conditions  by 
employing  a  single  component  scattering  object  that  exactly  satisfied 
the  (first-order)  Born  approximation.  The  propagation  of  determin¬ 
istic  artifacts  by  the  reconstruction  algorithms  under  less-than-ideal 
conditions  was  investigated  by  employing  a  two  component  scatter¬ 
ing  object  that  introduced  strong  and  multiple-scattering  effects  into 
the  measurement  data. 

A.  1  Single-Scattered  (Bom)  Data  and  Noise  Model  The  scat¬ 
tering  object  function,  shown  in  figure  4,  was  taken  to  be  a  lossless, 
uniform  cylindrical  disk  with  a  diameter  of  30  pixels  that  was 
convolved  with  a  symmetric  Gaussian  function  with  a  standard 
deviation  of  0.2  pixel.  The  Fourier  transform  of  the  object  function 
was  therefore  approximately  bandlimited  to  a  circular  disk  of  radius 
V2  v0  in  its  2D  Fourier  space.  It  was  assumed  that  the  scatterer  was 
weakly  scattering,  so  that  the  Bom  approximation  may  reasonably 
be  taken  to  hold.  Therefore,  these  numerical  simulations  were  de¬ 
signed  to  investigate  the  statistical  properties  of  the  reconstruction 
algorithms  rather  than  the  weak  scattering  model.  Minimal-scan  data 
sets  were  generated  by  using  the  FDP  theorem  to  calculate  simulated 
scattered  field  data  for  96  measurement  angles  <p  that  were  evenly 
spaced  between  0  and  37t/2.  At  each  measurement  angle,  128 
samples  were  calculated  with  a  sampling  increment  A£  =  1/2 v0, 
where  vQ  is  the  frequency  of  the  incident  plane  wave. 


The  integrals  in  the  curly  braces  of  the  first  and  second  terms  of 
equation  10  can  be  evaluated  (Metz  and  Pan,  1995)  yielding  the 
expressions  2  7 rjk  y  (■  vm)keik  eJk(  2  7rr  V  vfn  — ~vfj  and 

2 t r(-  1  )kjky( vm) " ke*keJk( 2 irrV vf„  v*  ),  respectively.  Using 
this  result  and  the  change  of  variables  va  =  Vv*  -  (which 
implies  vadvcl  =  v0/v'  vmdvm ),  we  arrive  at 


*  f  >£» 

a(r ,  0)  =  7 r  ^  /  [y kM’k(v„)  +  (-l)*y  kM,k(-v J] 

*— *  Vj-0 

X  ejkeJk(lTTvar)  va  dv(n  (11) 


where  vm  =  1  -  v2J 2v\.  Comparison  of  Eqns.  11  and  6b 

indicates  that,  for  0  <  va  ^  V2v0, 

PlM  =  \  [7 kM’k(v„)  +  (-1)V*A^("’ v»)].  (12) 


Figure  4.  The  original  scattering  object  function  was  formed  by 
convolving  a  uniform  circular  disk  with  a  diameter  of  30  pixels  with  a 
circularly  symmetric  Gaussian  function  with  a  standard  deviation  of 
0.2  pixels. 
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Table  I.  Error  values  of  the  reconstructed  images  shown  in  Figs.  9-11. 


Contrast 

MS-E-C  Error 

MS-FBPP  Error 

1.01 

12.03 

13.68 

1.05 

51.26 

53.70 

1.08 

115.25 

119.13 

V.  RESULTS 

A.  Single-Scattered  (Bom)  Data  Case.  We  first  used  the  MS- 
E-C  and  MS-FBPP  algorithms  to  reconstruct  the  scattering  object 
function  using  the  simulated  noiseless  minimal-scan  data  set.  The 
images  reconstructed  using  the  MS-E-C  and  MS-FBPP  algorithms 
are  displayed  in  figures  5a  and  5b,  respectively.  It  is  observed  that, 
in  the  absence  of  noise,  both  the  MS-E-C  and  MS-FBPP  algorithms 
can,  with  high  fidelity,  reconstruct  the  original  scattering  object 
function  from  the  minimal-scan  data  set. 

Using  one  of  the  noisy  minimal-scan  data  sets,  we  used  the 
MS-E-C  and  MS-FBPP  algorithms  to  reconstruct  the  scattering 
object  function.  The  noisy  images  reconstructed  using  the  MS-E-C 
and  MS-FBPP  algorithms  are  displayed  in  figures  6a  and  6b,  re¬ 
spectively.  The  image  reconstructed  using  the  MS-E-C  algorithm 
(Fig.  6a)  appears  less  affected  by  the  data  noise  and  more  closely 
resembles  the  original  object  than  does  the  image  reconstructed 
using  the  MS-FBPP  algorithm  (Fig.  6b).  The  local  image  variance, 
which  was  empirically  calculated  from  the  two  sets  of  250  noisy 
images  reconstructed  using  the  MS-E-C  and  MS-FBPP  algorithms, 
quantitatively  confirms  this  observation.  Figure  7  is  a  plot  of  the 
local  variance  obtained  from  the  MS-FBPP  reconstructed  images 
divided  by  the  local  variance  obtained  from  the  MS-E-C  recon¬ 
structed  images.  Clearly,  the  ratio  of  the  variances  is  everywhere 
greater  than  one,  and  near  the  corners  of  the  reconstructed  image  is 
as  great  as  ten.  This  quantitatively  demonstrates  that  the  MS-E-C 
reconstruction  algorithms  are  less  susceptible  to  the  effects  of  data 
noise  than  are  the  MS-FBPP  reconstruction  algorithms. 

B.  Multiple-Scattering  Case.  Using  the  MS-E-C  and  MS- 
FBPP  algorithms  we  reconstructed  the  two  component  scattering 
objects  shown  in  figures  9-11,  which  correspond  to  the  cases  where 
the  cylinders  had  refractive  index  values  of  n(r)  =  1.01,  1.05,  and 
1.08,  respectively.  In  each  case,  the  image  reconstructed  by  use  of 
the  MS-E-C  algorithm  (Figs.  9a-lla)  appears  to  contain  less  pro¬ 
nounced  artifacts  than  does  the  image  reconstructed  by  use  of  the 
MS-FBPP  algorithm  (Figs.  9b-llb).  This  observation  is  confirmed 
by  Table  I,  which  shows  that  for  each  value  of  n(r)  the  MS-E-C 
algorithm  produced  images  that  have  lower  error  values  than  the 
corresponding  images  produced  by  the  MS-FBPP  algorithm.  This 
quantitatively  demonstrates  that  the  MS-E-C  algorithms  are  less 
susceptible  to  multiple-scattering  effects  and  other  deterministic 
inconsistencies  than  are  the  MS-FBPP  algorithms.  However,  as  one 
would  expect,  the  performance  of  both  algorithms  dramatically 
deteriorates  as  the  refractive  index  values  increases  and  the  Born 
condition  (Chen  and  Stamnes,  1998)  is  severely  violated. 

VI.  DISCUSSION 

Previously  we  have  shown  (Pan  and  Anastasio,  1999)  that,  in  DT 
employing  the  2D  classical  scanning  geometry,  the  minimal-scan 
data  set  acquired  using  view  angles  only  in  [0,  3tt/2]  contains  all  of 
the  information  necessary  to  exactly  reconstruct  the  scattering  object 
function.  We  subsequently  developed  a  class  of  MS-FBPP  algo¬ 


rithms  that  were  capable  of  exactly  reconstructing  the  scattering 
object  function  from  the  minimal-scan  data  set. 

In  this  work,  we  have  developed  a  novel  class  of  reconstruction 
algorithms  for  the  minimal-scan  DT  reconstruction  problem.  These 
algorithms,  referred  to  as  MS-E-C  reconstruction  algorithms,  have 
distinct  advantages  over  the  MS-FBPP  reconstruction  algorithms. 
Because  the  FBPJ  algorithm  used  by  the  MS-E-C  algorithms  does 
not  involve  a  depth-dependent  filtering,  the  MS-E-C  algorithms  are 
more  computationally  efficient  than  are  the  MS-FBPP  algorithms. 
More  importantly,  we  have  quantitatively  demonstrated  that  the 
MS-E-C  algorithms  are  less  susceptible  to  data  noise,  modeling 
errors  due  to  the  violation  of  weak  scattering  conditions,  and  other 
finite  sampling  effects  than  are  the  MS-FBPP  algorithms.  This  result 
is  consistent  with  the  observation  that  the  MS-FBPP  algorithms 
involve  more  complicated  numerical  operations  than  do  the  MS-E-C 
algorithms,  which  may  amplify  the  propagation  of  noise  and  errors 
into  the  reconstructed  image.  Therefore,  the  use  of  a  MS-E-C 
algorithm  instead  of  a  MS-FBPP  algorithm  (using  the  same  weight¬ 
ing  function)  will  generally  result  in  a  reduction  of  the  reconstructed 
image  variance  and/or  a  reduction  of  the  image  artifacts. 

Recently,  non-linear  reconstruction  algorithms  that  incorporate 
higher-order  scattering  approximations  have  been  proposed  for  full- 
scan  DT  (Lu  and  Zhang,  1996;  Tsihrintzis  and  Devaney,  2000b; 
Tsihrintzis  and  Devaney,  2000a).  The  generalization  of  these  works 
to  case  of  minimal-scan  DT  is  an  important  task  that  is  currently 
under  way. 
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Abstract 

Diffraction  tomography  (DT)  is  an  inversion  technique  that  reconstructs  the  refractive  index  distribution  of  a  weakly  scatter¬ 
ing  object.  In  this  paper,  a  novel  reconstruction  algorithm  for  3D  diffraction  tomography  employing  spherical-wave  sources  is 
mathematically  developed  and  numerically  implemented.  Our  algorithm  is  numerically  robust  and  is  much  more  computationally 
efficient  than  the  conventional  filtered  backpropagation  algorithm.  Our  previously  developed  algorithm  for  DT  using  plane-wave 
sources  is  contained  as  a  special  case. 
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I.  Introduction 

In  diffraction  tomography  (DT),  a  weakly  scattering  object  is  interrogated  using  a  diffracting  wavefield, 
and  the  scattered  wavefield  around  the  object  is  measured  and  used  to  reconstruct  the  (low-pass  filtered) 
refractive  index  distribution  of  the  scattering  object.  The  principles  of  DT  have  been  extensively  utilized  for 
developing  optical  [1,2]  and  acoustic  [3]  tomographic  imaging  systems  for  biomedical  applications. 

It  is  widely  known  that  the  filtered  backpropagation  (FBPP)  and  direct  Fourier  (DF)  reconstruction  algo¬ 
rithms  for  three-dimensional  (3D)  DT  possess  certain  limitations  [4].  The  depth-dependent  filtering  (back- 
propagation)  in  the  3D  FBPP  algorithm  requires  a  large  number  of  2D  fast  Fourier  transforms  (FFTs)  to  be 
performed  for  processing  the  measured  data  at  each  measurement  view,  which  renders  the  3D  FBPP  algo¬ 
rithm  extremely  computationally  demanding.  Furthermore,  we  have  shown  that  (in  two-dimensional  (2D) 
DT)  the  FBPP  algorithm  may  considerably  amplify  data  noise  [5].  The  3D  DF  algorithms  require  the  use 
of  a  3D  interpolation  method  to  obtain  samples  on  a  3D  Cartesian  grid  in  the  Fourier  space  of  the  scattering 
object,  upon  which  a  3D  inverse  FFT  can  be  employed  to  reconstruct  the  scattering  object  function.  Because 
the  sample  density  in  the  3D  Fourier  space  obtained  from  the  measured  data  is  non-uniform,  sophisticated 
and  computationally  demanding  interpolation  strategies  are  generally  required  to  avoid  producing  significant 
interpolation  errors  that  would  degrade  the  accuracy  of  the  reconstructed  image. 

For  3D  DT  employing  plane-wave  sources,  we  have  recently  developed  a  new  class  of  reconstmction 
algorithms  that  circumvent  the  shortcomings  of  the  3D  FBPP  and  DF  algorithms  [4].  These  algorithms, 
referred  to  as  plane-wave  estimate-combination  (E-C)  reconstruction  algorithms,  effectively  reduce  the  3D 
DT  reconstruction  problem  to  a  series  of  2D  X-ray  reconstruction  problems,  and  thus  greatly  reduce  the 
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large  computational  load  required  by  conventional  3D  DT  reconstruction  algorithms.  Additionally,  these 
algorithms  do  not  require  an  explicit  3D  interpolation  in  the  Fourier  space  of  the  scattering  object. 

In  many  imaging  applications  [1,6],  it  may  be  useful  to  utilize  a  diverging  spherical-wave  rather  than  a 
plane-wave  to  interrogate  the  scattering  object.  Because  of  the  distinct  advantages  of  the  E-C  reconstruction 
algorithms  for  plane-wave  DT  [4],  it  is  important  to  generalize  them  to  DT  employing  spherical-wave  sources. 
In  this  work,  we  generalize  our  previously  developed  (plane-wave)  E-C  reconstruction  algorithms  to  DT 
employing  spherical-wave  sources  and  numerically  demonstrate  the  developed  algorithm  using  simulated 
data. 


II.  Background 

We  will  utilize  the  model  of  spherical- wave  DT  described  by  Devaney  in  [7],  the  scanning  geometry  of 
which  is  shown  in  Fig.  1.  The  case  of  3D  DT  utilizing  a  plane- wave  source  can  be  viewed  as  a  special  case 
of  spherical-wave  DT  and  will  be  discussed  below.  The  scattering  object  is  illuminated  by  monochromatic 
spherical -wave  source  located  at  the  position  rj  =  —  S  on  the  77-axis,  emitting  a  wavefield  of  the  form 


=  A0 


e<727ri/o|r— 57?| 

\r-sM 


=  Ao 


£2+z2+ (S-fD)2 

y^  +  ^-KS  +  D)2’ 
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where  A0  is  the  complex  amplitude,  k  —  2nu0  is  the  wavenumber,  /?  is  a  unit  vector  pointing  along  the 
positive  77-axis,  and  D  is  the  distance  of  the  detector  surface  from  the  center  of  rotation.  From  the  scattered 
wavefield  measured  in  the  £-z  plane  oriented  at  a  measurement  angle  </>  and  located  at  a  distance  77  =  D  from 
the  origin,  one  seeks  to  reconstruct  the  scattering  object  function  a{r).  The  underlying  physical  property  of 
the  scattering  object  that  is  mapped  in  DT  is  the  refractive  index  distribution  n(f),  which  is  related  to  the 
scattering  object  function  by  the  equation  a(r)  =  n2(r)  —  1. 

Let  u(£,  2,  <f>)  denote  the  measured  total  wavefield  in  the  £-z  plane  positioned  at  77  =  D,  as  shown  in  Fig. 
1.  The  scattered  data  is  given  by 


u,(£,  2, 4)  =  u(£,  z,  <{> )  -  u0(£,  4>), 

which  can  be  treated  as  a  measurable  data  function  because  u(£,  z,  <f> )  can  be  measured  and  because  uo  (£,  2, 
is  assumed  known.  We  can  introduce  a  modified  data  function 
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and 
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Using  paraxial1  and  weak  scattering  approximations,  Devaney  derived  [7]  the  spherical-wave  FDP  theorem 
that  is  given  by 


O00 

a(r)< 

-00 


rjWept+ff+V-nWtf  if  v2m  +  u2  <  X2U2 

-00  J  — 00 

^he  paraxial  approximation  requires  that  both  S  and  D  are  much  larger  than  the  size  of  the  scattering  object. 
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=  0  iivl  +  vl>X2vl  (6) 

Equation  6  states  that  the  modified  data  function,  M (vm,  i /*,</>),  is  equal  to  a  semi-ellipsoidal  slice,  oriented 
at  angle  <j>,  through  the  3D  Fourier  transform  of  the  object  function  a(r).  In  this  work  we  assume  the  validity 
of  the  first-order  Bom  approximation;  the  spherical-wave  FDP  theorem  can  also  be  derived  based  upon  the 
Rytov  approximation  by  appropriately  redefining  Eqn.  3  [7]. 

Equation  6,  coupled  with  a  3D  interpolation  method,  can  be  used  to  implement  a  DF  reconstruction  algo¬ 
rithm.  The  3D  FBPP  algorithm  for  plane-wave  DT  [8]  can  readily  be  extended  to  3D  spherical-wave  DT  [6] 
and  is  given  by  [6, 9] 


(7) 


where 


v,  =  -  x>K2  +  K2)  ~  •*)■ 

X  X 


III.  E-C  Reconstruction  Algorithm 

In  deriving  the  E-C  reconstructions  algorithms  for  spherical-wave  DT,  we  will  modify  the  procedure  em¬ 
ployed  for  the  derivation  of  the  plane- wave  DT  algorithms  described  in  [4].  Let  a(r,  6,  z)  denote  a  3D  scat¬ 
tering  object  function  in  the  cylindrical  coordinate  system.  The  X-ray  transform,  p  (£,  z,  of  a(r,  9,  z)  is 
defined  as 

roc 

P  (£,  2,  <t>a)  =  /  a{r,  0,  z)dr],  (8) 

where  <f>a  is  the  projection  angle,  £  =  rcos(</»o  —  0)  and  77  =  —  rsin(</>0  —  9).  Equation  8  states  that  p(£,  z,  <f>o)  is 
the  geometrical  projection  of  a(r,  9 ,  z)  onto  the  £— z  plane  oriented  at  angle  (j)0  about  the  z  axis.  Consequently, 
p(£,  z,  <f>o)  can  be  interpreted  as  a  stack  of  2D  Radon  transforms  of  a(r,  9 ,  z)  along  the  z  axis.  The  combination 
of  a  2D  Fourier  transform  with  respect  to  £  and  z  and  a  ID  Fourier  series  expansion  with  respect  to  </>0  of 
p(£,  z,  <t>0)  is  given  by 

1  /*2tt  /*  roo 

pk(^,^)=f  //  (9) 

Z  7T  J^z——oo 


where  the  integer  k  is  the  angular  frequency  index  with  respect  to  </>0,  and  va  and  vz  are  the  continuous  spatial 
frequencies  conjugate  to  £  and  z,  respectively.  As  a  matter  of  convenience,  we  refer  to  Pfc(VQ,  vz)  as  the  “3D 
Fourier  transform”  of  p  (£,  z,  4>q).  Substituting  Eqn.  8  into  Eqn.  9  and  carrying  out  the  integral  over  (f> 0  yields 
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Again  for  convenience,  we  refer  to  vz)  =  ^  /027r  M(um,  uz,  (f>)  d(f>  as  the  “3D  Fourier  trans¬ 

form”  of  the  modified  data  function.  Using  Eqn.  6  and  noting  that  £  =  r  cos (<f>  —  9)  and  tj  =  —r  sin(<^  —  9) 
one  can  obtain 


{ if 
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g  27 Tv^r  sin( 


r^rcos(^-fl))  ^  j  r  drdQ, 


(ID 


where  v ^  +  <  X2rf-  Carrying  out  the  integral  [10]  in  the  curly  brackets  in  Eqn.  11,  one  can  re-express 

Mk(um,uz)  as 
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Mk(um,u.)  =Hy[-£+2M‘r  .-i**  <b 

\  Vf  —  V2  JZ=- oo 
Y  m  /i 

noo  _  / - - - 

a(r,  0, 2)  e_jW  3k(2nr^/i^2  -  v*)  r  dr  dd 

=0  VP 


where  v\  +  i >\<  x2v$. 

Comparison  of  Eqns.  10  and  12  yields  that  for  +  v2  <  x2^o 


P k(Va,  V'z)  =  [7 Wmi  V'z)]k  Mk(^m,  V,)> 


(12) 


(13) 


where 

and 


"l  =  V'm 


2~»l 


7K>»  v'z) 


.III1  2  —  ifi- 
Y  " m  11 

VL+Vu 


(14) 

(15) 


The  explicit  forms  of  this  general  relationship  in  Eqn.  13  are  determined  completely  by  the  forms  of  the 
relationship  between  ua  and  vm  implied  in  Eqn.  14.  Equation  14  indicates  that  the  condition  "I  +  "I  <  X2< 
is  equivalent  to  the  condition  v2  +  v'2  <  1  +  1/x2).  In  order  to  determine  vm  as  a  function  of  vz  and  ua, 

we  must  solve  for  the  roots  of  the  fourth-order  equation 


Ci +  Civ’m  +  C3  —  0, 


(16) 


where  the  coefficients  Ci  are  given  by 


Ci  =  (1  -  x2)2, 


vt 


C2  =  2(1  —  X2)(2^o  ~  ua  ~  ~o)  +  4i'oX2> 


and 


X" 


C3  =  (2u2  -  n2  -  Vl)2  -  4 v2M  -  4. 

A  A 


The  four  roots  of  Eqn.  16  are  given  by 


"'ml  = 


m2 


-c2  +  v fc2  -  4C1C'3 

2  Ci 


1/2 


(17) 


and 


Um3  ~  ~vmi  ~ 


-c2  -  Jcl  -  4CXC3 

2C\ 


1/2 


(18) 


Let  vmi  —  i  =  1,2, 3,  and  4.  For  v\  +  v'2  <  Vq(1  +  l/x2)>  the  two  roots  v'mZ  and  u'mA  are  complex¬ 

valued  and  therefore  not  physically  meaningful.  In  the  plane-wave  case  (x  =  1),  there  are  only  two  roots  that 
are  given  by  [4] 

Vml  =  -Vm2  =  I'm  =  { (l£  +  Itj)  [l  -  ^2^1  ~  '  >  (19) 

where  n2  +  v2  <  2nz. 
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Therefore,  for  a  given  pair  of  va  and  v'z  satisfying  0  <  + 1/2  <  u^(l  +  1/x2),  Pk(va,v'z)  can  be  obtained 

from  Mfc (i/m,  uz)  at  vmX  =  um  and  uz  as 

=  frK.,  ^)]fcMfc(l/m,  (20) 

and  also  from  Mk{vm,  vz)  at  i/m2  =  —vm  and  vz  as 


P h(Va,  v'z)  =  h(um2,  Vz)}kMk(vm2,  Vz) 

=  (-^)kh(v'm,v’z)]~kMk(-um,uz).  (21) 

The  fact  that  there  are  two  distinct  ways  to  obtain  Pk(va,  K)  from  the  measured  data  can  be  explained  by  the 
existence  of  a  double  coverage  of  the  scattering  object  Fourier  space  [7],  At  a  given  measurement  angle  0, 
the  spherical-wave  FDP  theorem  relates  the  2D  Fourier  transform  of  the  modified  measured  data  to  the  3D 
Fourier  transform  of  a(r)  evaluated  over  a  semi-ellipsoid  oriented  at  angle  0.  As  the  measurement  angle  0 
is  varied  from  0  to  2tt,  two  overlapping  coverages  of  the  Fourier  space  are  generated,  with  each  coverage 
producing  one  of  the  relationships  described  by  Eqns.  20  or  21. 

Equations  20  and  21  yield  two  identical  values  of  Pfc(z^a,  v'z)  when  the  measured  data  are  consistent.  How¬ 
ever,  when  the  measured  data  contain  noise,  Eqns.  20  and  21  will  generally  produce  different  values  of 
Pfc(^a)  v'z).  These  two  values  can  be  combined  to  obtain  a  final  estimate  of  Pk  (i/a,  v'z)  that  has  a  reduced  noise 
level  as 

piW)  K,  v'z)  =  uk(ym,  vz)[ 7(z4,  v'z)}kMk{vm,  vz) 

+  (1  -  vz))(-\)k{^{v'm,  v'z)]-kMk(-i /m>  uz),  (22) 

where  u)k{vm,  vz)  is  a  generally  complex-valued  combination  coefficient.  The  superscript  “w”  indicates  that 
P v'z)  is  obtained  by  use  of  a  combination  coefficient  uk{ Vm,  vz).  If  Mfc(i/m,  vz)  and  Mfc(— i/m,  vz)  are 
interpreted  as  a  random  variables,  then  for  a  given  uk{vm,  vz),  Eqn.  22  can  be  interpreted  as  an  estimation 
method  for  obtaining  the  ideal  sinogram.  Because  oJk(vm,  vz)  may  be  any  complex- valued  function  of  um,  vz 
and  k,  Eqn.  22,  in  effect,  represents  an  infinite  class  of  estimation  methods.  An  estimate  of  p(£,  z,  (f>0)  can  be 
obtained  by  taking  the  3D  inverse  Fourier  transform  of  Pj^  (ua,  uz).  For  a  fixed  value  of  z,  the  filtered  back- 
projection  (FBP)  algorithm  of  X-ray  CT  can  be  employed  to  reconstruct  the  corresponding  transverse  slice 
of  a(r,  9,  z )  from  p(£,  z,  <fi0).  We  refer  to  the  combination  of  Eqn.  22  to  estimate  Pj^  (um,  uz)  coupled  with 
the  2D  FBP  algorithm  to  reconstruct  transverse  slices  of  a(r,  9,  z)  as  a  spherical-wave  E-C  reconstruction 
algorithm  for  3D  DT. 


IV.  Numerical  Results 

We  performed  numerical  simulations  to  demonstrate  the  spherical-wave  E-C  reconstruction  algorithm. 
We  considered  a  mathematical  phantom  comprised  of  two  different  (uniform)  spheres  whose  3D  Fourier 
transforms  were  approximately  bandlimited  to  a  sphere  of  radius  \/2u0.  Our  intention  was  not  to  test  the 
validity  of  the  weak  scattering  (Bom)  condition  (see  the  Discussion  Section),  but  rather  to  demonstrate  that 
the  spherical-wave  E-C  reconstruction  algorithm  can  accurately  reconstruct  the  scattering  object  function 
from  weakly  scattered  data.  We  therefore  employed  Eqns.  3  and  6,  along  with  the  analytic  expression  for 
the  3D  Fourier  transform  of  a  the  spheres,  to  calculate  noiseless  samples  of  us(£,  z,  <j>)  over  a  128  x  128 
detector  array  at  128  view  angles  that  were  evenly  spaced  over  360°.  In  generating  the  simulation  data, 
we  assumed  a  scanning  geometry  with  S=100  (arbitrary  units)  and  D=104.0816  that,  according  to  Eqn.  4, 
yields  x  =  0.7.  In  order  to  simulate  the  stochastic  nature  of  noisy  scattered  data,  we  created  a  second 
data  set,  us(£,  z,  f),  where  the  measured  scattered  data  were  treated  as  samples  of  an  uncorrelated  bivariate 
Gaussian  stochastic  process.  When  generating  us(£,  z,  0),  the  mean  and  variance  parameters  describing  the 
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real  (imaginary)  component  of  the  stochastic  process  were  set  equal  to  the  real  (imaginary)  values  of  the 
noiseless  data  us(f,  z,  <f>).  Therefore,  at  a  given  position  (£,  z,  <f>)  in  the  data  space,  the  magnitude  of  noise 
contained  in  the  real  and  imaginary  components  of  us(£,  z,  <f>)  was  proportional  to  the  magnitude  of  the  real 
and  imaginary  components  of  us  (£,  z,  <fi),  respectively. 

We  reconstructed  three  transverse  slices  of  the  scattering  object  by  use  of  the  spherical-wave  E-C  algorithm 
specified  by  tok{vm,  vz)  =  |  in  Eqn.  22.  It  can  readily  be  shown  that  when  the  data  noise  is  uncorrelated, 

u  =  |  is  a  statistically  optimal  choice  in  the  sense  that  it  minimizes  the  variance  of  the  estimate  Pj^(t'a), 
which  in  turn,  results  in  a  minimization  of  the  global  image  variance2  [4].  For  more  general  noise  models,  it  is, 
in  principle,  possible  to  derive  other  optimal  forms  for  uk(vm,  vz).  The  true  images  of  the  chosen  transverse 
slices  are  shown  in  Fig.  3-(a).  The  images  reconstructed  from  the  noiseless  data,  shown  in  Fig.  3-(b),  do 
not  contain  any  artifacts  and  accurately  represent  the  corresponding  true  slices  shown  in  Fig.  3-(a).  The 
same  images  reconstructed  from  the  noisy  data  set  are  shown  in  Fig.  3-(c).  Although  noisy  in  appearance, 
the  images  are  not  distorted  which  confirms  that  the  spherical-wave  E-C  algorithm  are  numerically  stable. 
Out  of  curiosity,  we  also  reconstructed  the  same  transverse  slices  by  use  of  our  previously  developed  plane- 
wave  E-C  reconstruction  algorithm  (that  assumes  the  measurement  geometry  corresponds  to  x  =  1)  and  the 
noiseless  data  set.  The  reconstructed  images,  shown  in  Fig.  4,  clearly  contain  artifacts  and  distortions.  This 
numerically  demonstrates  the  importance  of  properly  accounting  for  the  wavefield  curvature  in  3D  DT. 

V.  Discussion 

Previously,  we  developed  a  novel  class  of  reconstruction  algorithms  for  3D  DT  using  plane-wave  sources. 
These  algorithms,  referred  to  as  plane-wave  E-C  reconstruction  algorithms,  had  a  significant  computational 
advantage  over  the  conventional  3D  FBPP  algorithm,  and  unlike  the  3D  DF  method,  did  not  require  an  explicit 
3D  interpolation  in  the  Fourier  space  of  the  scattering  object. 

The  use  of  a  plane-wave  source  may  not  be  feasible  in  many  experimental  situations,  and  it  may  be  more 
convenient  to  interrogate  the  scattering  object  using  a  diverging  spherical  wave  that  is  produced  by  a  point 
source.  In  this  work,  we  have  developed  a  class  of  spherical- wave  E-C  reconstruction  algorithms  for  DT 
using  spherical-wave  sources  and  measurement  geometries  that  satisfy  the  paraxial  approximation  [7].  The 
spherical-wave  E-C  reconstruction  algorithms  can  be  viewed  as  generalizations  of  the  plane-wave  E-C  re¬ 
construction  algorithms,  and  reduce  to  the  plane- wave  E-C  algorithms  in  the  special  case  x  =  1.  The 
spherical-wave  E-C  reconstruction  algorithms  possess  the  same  Advantages  as  their  plane-wave  counterparts. 
For  example,  to  reconstruct  an  Nz  image  volume  the  spherical-wave  E-C  algorithm  requires  ~  jV3logiV  nu¬ 
merical  operations  while  the  spherical-wave  FBPP  algorithm  described  by  Eqn.  7  would  require  ~  N4\ogN 
numerical  operations.  Unlike  DF  methods,  the  spherical-wave  E-C  algorithms  do  not  require  an  explicit  3D 
interpolation  in  the  Fourier  space  of  the  scattering  object. 

The  spherical-wave  E-C  reconstruction  algorithms  have  been  developed  using  the  first-order  Bom  (or  Ry- 
tov)  weak  scattering  approximation.  In  certain  applications,  the  weak  scattering  approximation  may  not  be 
valid  and  the  reconstructed  image  may  contain  artifacts.  However,  the  spherical-wave  E-C  algorithms  provide 
a  natural  framework  for  the  incorporation  of  higher-order  scattering  perturbation  approximations  [11—13]  into 
the  algorithms. 
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Figure  Captions 

Fig.  1.  The  classical  scanning  geometry  of  spherical-wave  3D  DT.  A  spherical-wave  is  incident  at  angle  0, 
and  the  scattered  field  is  measured  in  the  £-z  plane  positioned  at  77  =  D.  The  measurement  angle  <p  is  varied 
from  0  to  27 r. 

Fig.  2.  A  plot  of  vm  as  a  function  of  vz  and  v0  as  described  by  Eqn.  17.  The  plot  was  generated  using  x  —  0.7 

and  i/0  =  7r. 


Fig.  3.  (a)  True  slices  through  the  numerical  phantom.  Transverse  slices  reconstructed  from  (b)  noiseless  and 
(c)  noisy  data  using  a  spherical-wave  E-C  reconstruction  algorithm. 

Fig.  4.  Transverse  slices  reconstructed  from  the  spherical- wave  data  function  (x  =  0.7)  using  a  plane- wave 
DT  reconstruction  algorithm.  As  expected,  the  images  contain  distortions. 
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Fig.  4. 
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Abstract 

In  reflectivity  tomography,  conventional  reconstruction  approaches  require  that  mea¬ 
surements  be  acquired  at  view  angles  that  span  a  full  angular  range  of  27r.  It  is  often, 
however,  advantageous  to  reduce  the  angular  range  over  which  measurements  are  acquired, 
which  can,  for  example,  minimize  artifacts  due  to  movements  of  the  imaged  object.  More¬ 
over,  in  certain  situations  it  may  not  be  experimentally  possible  to  collect  data  over  a  27r 
angular  range.  In  this  work,  we  investigate  the  problem  of  reconstructing  images  from 
reduced-scan  data  in  reflectivity  tomography.  By  exploiting  symmetries  in  the  data  func¬ 
tion  of  reflectivity  tomography,  we  heuristically  demonstrate  that  an  image  function  can  be 
uniquely  specified  by  reduced-scan  data  that  correspond  to  measurements  taken  over  an 
angular  interval  (possibly  disjoint)  that  spans  at  least  7r  radians.  We  also  identify  sufficient 
conditions  that  permit  for  a  stable  reconstruction  of  image  boundaries  from  reduced-scan 
data.  Numerical  results  in  computer-simulation  studies  indicate  that  images  can  be  recon¬ 
structed  accurately  from  reduced-scan  data. 

EDICS  Category:  3-TOMO  (Computer  tomography) 


1  Introduction 

Reflectivity  tomography  has  been  applied  to  numerous  biomedical  [1,2]  and  non-destructive 
testing  [3]  imaging  problems.  In  two-dimensional  (2D)  reflectivity  tomography  [4,  5],  a 
weakly  reflecting  object  that  is  immersed  in  an  acoustically  homogeneous  background 
medium  is  illuminated  with  infinitesimally  short  ultrasonic  pulses,  and  the  concomitant 
reflected  signals  are  measured  as  functions  of  time  at  each  of  the  multiple  source  locations. 
The  task  in  reflectivity  tomography  is  to  reconstruct  from  such  measured  data  a  function 
describing  the  reflectivity  distribution  within  the  scattering  object1. 

In  reflectivity  tomography,  an  imaging  configuration  that  acquires  data  over  the  com¬ 
plete  angular  interval  of  27r  is  referred  to  as  a  full-scan  configuration.  Norton  derived  [4] 
an  analytic  algorithm  for  reconstructing  the  object  function  from  data  acquired  with  a 
full-scan  configuration.  In  tomographic  imaging,  however,  it  is  often  desirable  to  minimize 
the  angular  range  over  which  data  are  acquired.  When  ionizing  radiation  is  employed, 
this  can  reduce  the  radiation  dose  that  is  delivered  to  the  subject  during  the  imaging  pro¬ 
cess.  An  advantage  of  reducing  the  scanning  angle,  which  is  also  relevant  to  biomedical 
imaging  applications  employing  non-ionizing  radiation,  is  that  the  time  needed  to  collect 
data  is  reduced,  thus  diminishing  image  artifacts  and  distortions  created  by  movements  of 

1For  simplicity,  we  refer  to  the  object’s  reflectivity  function  as  the  object  function. 
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the  object  during  the  imaging  process.  Furthermore,  in  certain  situations  it  may  not  be 
experimentally  possible  to  collect  data  over  a  full  angular  range  of  2ir. 

An  imaging  configuration  that  acquires  data  over  an  angular  interval  <f>,  where  7r  < 
$  <  27t,  is  referred  to  as  a  reduced-scan  configuration.  It  has  been  shown  that  one  can 
accurately  reconstruct  images  from  data  acquired  with  reduced-scan  configurations  in  to¬ 
mographic  imaging  modalities  including  fan-beam  computed  tomography  (CT)  [6]  and 
single-photon  emission  computed  tomography  (SPECT)  [7,8].  However,  it  remains  un¬ 
clear  whether  an  accurate  object  function  can  be  reconstructed  from  reduced-scan  data  in 
reflectivity  tomography. 

In  this  work,  we  investigate  the  problem  of  image  reconstruction  from  reduced-scan 
data  in  reflectivity  tomography.  A  potato-peeler  perspective2  is  proposed  for  identifying 
data  symmetries  in  reflectivity  tomography.  Using  the  identified  data  symmetries,  we 
heuristically  demonstrate  that  data  acquired  with  reduced-scan  configurations  in  reflectivity 
tomography  are,  in  principle,  sufficient  to  specify  the  object  function.  We  also  apply  results 
from  microlocal  analysis  [10]  to  investigate  the  stability  of  reconstructing  image  boundaries 
from  reduced-scan  data  in  reflectivity  tomography. 

Norton’s  algorithm  [4],  which  can  be  applied  only  to  full-scan  data,  cannot  accurately 
reconstruct  images  from  reduced-scan  data.  Instead,  we  propose  to  use  the  expectation 
maximization  (EM)  algorithm  [11]  for  reconstruction  of  the  object  function  from  reduced- 
scan  data.  Using  computer-simulation  studies,  we  evaluate  the  numerical  properties  of 
images  that  are  reconstructed  from  reduced-scan  data.  The  results  of  these  studies  indicate 
that,  under  practically  relevant  conditions,  images  can  be  reconstructed  accurately  from 
reduced-scan  data.  In  particular,  for  reflecting  objects  and  scanning  configurations  that 
may  arise  in  realistic  experiments,  numerically  accurate  images  can  be  reconstructed  from 
reduced-scan  data  acquired  over,  e.g.,  an  angular  interval  of  only  7r. 

2  Data  Function  in  Reflectivity  Tomography 

In  this  section,  we  briefly  discuss  the  data  function  in  reflectivity  tomography  and  summa¬ 
rize  Norton’s  algorithm  [4]  for  image  reconstruction  from  full-scan  data.  Let  /(r,  9)  denote  a 
two-dimensional  (2D),  bounded,  and  non-negative  object  function  with  compact  support  on 
the  disk  of  radius  Rf  centered  at  the  origin,  for  some  fixed  real  number  Rf  >  0.  As  shown 

2The  potato-peeler  perspective  is  conceptually  similar  to  the  so-called  “layer-stripping”  method  for 
solving  the  Helmholtz  equation,  e.g.,  in  impedance  imaging  [9],  in  the  sense  that  the  object  function 
is  determined  layer  by  layer.  However,  the  mathematical  formulation  of  the  potato-peeler  perspective 
differs  entirely  from  that  of  the  existing  “layer-stripping”  method  [9].  We  are  currently  developing  a 
mathematical  formulation  for  the  potato-peeler  perspective  for  reflectivity  tomography  and  will  report 
such  results  elsewhere. 
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in  Fig.  1,  the  object  function  is  embedded  in  an  acoustically  homogeneous  background 
medium,  and  an  omni-directional  point  source-receiver  is  located  at  the  polar  coordinates 
(R,  </>  +  7r)  on  a  circular  source-receiver  trajectory  curve  C  that  encloses  the  object  function. 
At  time  t  =  0,  an  infinitesimally  short  pulse  of  sound  is  emitted  from  the  point  source,  and 
the  wavefield  that  is  backscattered  (i.e.,  reflected)  from  the  reflecting  object  is  measured 
at  the  source  location  as  a  function  of  time. 

The  measured  data  can,  under  certain  conditions  [5],  be  interpreted  as  a  set  of  line 
integrals  of  f(r,9)  defined  over  all  circular  arcs  (with  varying  radii)  that  are  centered  at 
a  source-receiver  location  (R,  <f>  +  7r),  where  Rf  <  R.  If  this  transmit/receive  process  is 
repeated  for  all  points  on  C,  multiple  sets  of  these  line  integrals  are  obtained.  These  line 
integrals  can  be  expressed  mathematically  as 

p(^i</*)  =  J  d9  J  f  drr  f(r,  6)  9  +  r2  —  2i?r  cos((/>  —  0)  —  R  —  £),  (1) 

where  £  G  [~R,R],  and  the  Dirac  5  function  restricts  the  integration  over  a  circular  arc 
with  a  radius  R  +  £  >  0  and  centered  at  the  point  (R,  <f>  +  7r)  on  the  curve  C.  The 
reconstruction  problem  in  reflectivity  tomography  is  to  determine  the  object  function  / (r,  6) 
from  knowledge  of  the  data  function  g(£,  (f>). 

In  full-scan  reflectivity  tomography,  because  #(£,  4>)  is  measured  for  all  <f>  G  [0,  27t),  one 
can  calculate  its  Fourier  series  expansion  coefficient  as 

(2) 

The  object  function  can  be  expressed  as  a  Fourier  series, 

OO 

/M)=  £  (3) 

n— — oo 


where  fn(r)  =  ^  /02ir  d9  f(r,  6)  e~jn6.  Norton  [4]  showed  that  fn(r)  (i.e.,  the  object  function) 
can  be  determined  from  knowledge  of  gn(£)  (i.e.,  the  full-scan  data  </(£,  </>))  as 


f  (r)  -  f°°dzz  3"{tz) 
M  '  Jo  J„(Rz) 


i 


R  d^9n^J°^ 


R 


2np 


(4) 


where  p  =  R  +  £  and  Jn(-)  is  the  nth-order  Bessel  function  of  the  first  kind.  Therefore, 
Eqns.  2,  3,  and  4  provide  a  recipe  for  reconstruction  of  the  object  function  from  full-scan 
data. 


3  Data  Symmetry  and  Reduction  of  Scanning  Angle 

Data  symmetries  have  been  exploited  to  demonstrate  whether  accurate  object  functions 
can  be  reconstructed  from  reduced-scan  data  in  other  tomographic  imaging  modalities 
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(e.g.,  fan-beam  CT  [6]  and  SPECT  [7,8]).  Below,  we  use  the  Radon  transform  of  an  object 
function  as  a  data  function  to  illustrate  a  data  symmetry.  The  Radon  transform  of  an 
object  function  f(r,9)  is  defined  as 

r2ir  rRf 

p(£A)=  \  dd  drr  f(r,  9)  6(£  —  rcos((f>  —  9)),  (5) 

Jo  Jo 

where  £  €  [—/?/,  Rj],  and  <j>  is  the  projection  angle.  A  symmetry  of  the  Radon  transform 
exists  because  its  values  at  conjugate  views  are  identical,  i.e., 

p{(,</>)  =P(~Z,  <(>  +  *)■  (6) 

Such  a  symmetry  suggests  that  the  full-scan  Radon  transform  contains  redundant  infor¬ 
mation.  Knowledge  of  the  Radon  transform  acquired  over  [0,  n)  can  be  used  for  recon¬ 
structing  an  image  exactly  because  knowledge  of  the  Radon  transform  corresponding  to 
the  un-measured  angular  intervals  is  fully  specified  by  Eqn.  6.  Inspection  of  Eqn.  1  in¬ 
dicates  that  the  data  function  g  (£,</>)  in  reflectivity  tomography  possesses  no  symmetry 
analogous  to  that  of  the  Radon  transform  in  Eqn.  6.  However,  we  reveal  below  that  the 
data  function  g(£,  4>)  does  admit  symmetries  and  that  such  symmetries  can  be  exploited  to 
identify  redundant  information  in  full-scan  data  in  reflectivity  tomography. 

3.1  Potato-Peeler  Perspective  for  Reflectivity  Tomography 

Although  the  data  function  in  Eqn.  1  does  not  admit  a  simple  symmetry  analogous  to 
that  of  the  Radon  transform  in  Eqn.  6,  the  fact  that  the  imaging  transform  in  Eqn.  1  is 
linear  and  that  f(r,  9)  has  a  compact  support  can  facilitate  the  explicit  identification  of 
data  symmetries  in  reflectivity  tomography.  We  generalize  the  potato-peeler  perspective, 
which  was  developed  previously  for  analyzing  data  symmetry  in  SPECT  [12],  to  investigate 
data  symmetries  in  reflectivity  tomography.  It  should  be  emphasized  that  this  perspective, 
to  be  presented  below,  is  not  intended  as  a  mathematically  stable  and  computationally 
practical  reconstruction  algorithm.  Instead,  it  is  used  only  for  heuristically  revealing  data 
symmetries  in  reflectivity  tomography. 

Without  loss  of  generality,  as  shown  in  Fig.  2,  we  assume  that  the  non-zero  support  of 
/(r,  9)  is  a  disk  of  radius  Rf  centered  at  the  origin,  i.e.,  f(r,  9)  ^  0  for  r  <  Rj  and  /(r,  9)  =  0 
for  r  >  Rf.  (The  procedures  and  discussion  to  be  presented  below  can,  however,  readily  be 
modified  to  address  any  compactly  supported  object  functions.)  We  use  L(£,  4>)  to  denote 
the  circular  integrating  line  in  Eqn.  1. 

At  a  given  view  angle  <fio,  one  can  identify  a  value  £mai  >  0  such  that  g(£max,(f>o)  /  0 
and  <?(£,</> o)  =  0  for  £  >  £mai.  The  assumption  on  the  support  compactness  of  f(r,9) 
indicates  that 

£mai  =  dij.  (7) 
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Therefore,  one  can  determine  Rf  from  knowledge  of  £max.  Conceptually,  we  treat  /(r,  9)  as 
a  collection  of  weighted  2D  Dirac  delta  functions.3  Because  the  integrating  line  L(£max,  0O) 
intersects  fir, 9)  only  at  the  single  point  (Rf,4> o)  (i  e.,  point  B  in  Fig.  2a),  Eqn.  1  reduces 
to 

9ifmax ,  <f> o)  =  f(Rf  i  0o)-  (8) 

On  the  other  hand,  at  the  conjugate  view  angle  0 o  +  7r,  one  can  identify  a  same  value 
Uax  >  o  such  that  g(-imax,  4>0+n)  ^  0  and  g{-£,  0o+tt)  =  0  for  £  >  fmax •  In  this  situation, 
the  integrating  line  L(—£max,(j)  o  +  vr)  also  intersects  only  at  a  single  point  (Rf,  <^o)  (i-e., 
point  B  in  Fig.  2a)  on  the  outermost  edge  of  the  support.  Again,  as  a  result,  Eqn.  1 
reduces  to 

g{  f,maxi  00  T  71”)  /(-^/>0 o)-  (9) 

Therefore,  the  object  function  f(r,0 )  at  the  outermost  point  (Rf,  0O)  (i.e.,  point  B  in  Fig. 
2a)  can  be  specified  completely  by  either  one  of  the  two  measurements  <?(£max,  0o)  and 
g(  £max j  00  "b  7T.)- 

The  component  of  f(r,9)  corresponding  to  the  point  (Rf,(j) 0)  produces  (via.  Eqn.  1)  a 
component  of  the  data  function  given  by 

gB(fiA)  =  9(Uax,<t>zmA),  (10) 

where  /(£,  0)  is  an  indicator  function  defined  as  /(£,  0)  =  1  for  £  =  [R2j +R?— 2RfRcos (0o — 
0)]  2  —  R  and  /(£,  0)  =  0  otherwise.  As  shown  in  Fig.  2b,  the  function  /(£,  0)  specifies 
the  loci  of  points  in  the  data  space  that  receive  a  contribution  when  f(r,  9)  in  Eqn.  1  is 
replaced  by  f(Rj,  0o)l?O'~R/)ftg~M-  From  g(f,  0)  and  #b(£,  0),  one  can  obtain  a  new  data 
function  as 

5(f,0)  =  0(f>0)-  0b(£,0),  (11) 

which  contains  no  contributions  from  point  B  in  the  object  function.  In  a  similar  manner, 
all  points  (e.g.,  points  A  and  C)  on  the  outermost  edge  of  f(r,9 )  can  be  determined  and 
“peeled  away”,  and  their  contributions  to  the  data  function  ^(^,0)  can  be  removed,  each 
time  forming  a  new  data  function.  The  next  set  of  points  at  the  interior  of  the  object 
function  will  be  exposed.  By  the  same  procedure,  the  next  layer  of  the  object  function  can 
be  determined  and  “peeled  away” ,  and  their  contributions  to  the  new  data  function  can 
be  removed.  By  repeating  this  procedure,  one  can  determine  f(r,  9)  by  completely  peeling, 
layer  by  layer,  the  object  function  away. 

3i.e„  /(r,  6)  =  J f  iff  dr'r'  f(r',0')  Sir=!^i!=Sl, 
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3.2  Reduction  of  Scanning  Angle 

From  Eqns.  8  and  9,  one  can  readily  identify  a  symmetry  of  the  data  function  as 

9{£maxi  4*o)  di.  4>maxi  4*0  “F  (12) 

Furthermore,  using  the  potato-peeler  perspective  described  above,  one  can  identify  a  sym¬ 
metry  of  g(£,  4>)  as 

g(U4>)  =  g(-U4>  +  i r),  (13) 

where  g(£m,  4>)  ±  0  and  4>  +  tt)  /  0  for  0  <  £m  <  £ max ,  and  g(£,  4>)  =  0  for  |£|  >  £m. 

The  symmetry  condition  in  Eqn.  13  indicates  that  any  point  of  the  object  function  can 
be  specified  by  the  values  of  a  (processed)  data  function  <?(£,  4>)  at  conjugate  views.  As  a 
result,  the  potato-peeler  perspective  suggests  that  the  object  function  may  be  specified  by 
knowledge  of  the  measured  data  g(^,4*)  available  over  </>  G  [0,7r).  We  refer  to  a  scanning 
configuration  that  acquires  g(£,  4>)  over  4>  €  [0, 7r)  as  the  short-scan  configuration.  More 
generally,  the  potato-peeler  perspective  indicates  that  an  object  function  may  be  specified 
by  knowledge  of  data  g(£  4 *)  measured  over  </>  G  4?,  where  $  is  any  proper  subset  of  [0, 2n) 
such  that  4>'  G  $  OR  </>'  +  7r  G  4>,  for  all  4>'  G  [0, 7r).  We  refer  to  a  scanning  configuration 
that  acquires  data  at  </>  G  <f>  as  a  it -scheme  configuration  [8].  Obviously,  the  short-scan 
configuration  can  be  viewed  as  a  special  case  of  the  7r-scheme  configurations.  The  data 
symmetries  above  indicate  that  full-scan  data  in  reflectivity  tomography  contain  redundant 
information  that  can  be  exploited  for  reducing  the  angular  range  over  which  measurements 
need  to  be  acquired. 

By  use  of  the  potato-peeler  perspective,  it  was  revealed  heuristically  above  that  an  ob¬ 
ject  function  f(r,9)  in  reflectivity  tomography  can  be  specified  uniquely  by  knowledge  of 
short-scan  or  7r-scheme  data.  This  is  consistent  with  the  result  of  an  analysis  of  a  family  of 
Radon  transforms  over  circles  [13].  One  can  use  principles  from  microlocal  analysis  [10,14], 
as  described  in  Appendix  A,  to  derive  theoretically  sufficient  conditions  for  stable  recon¬ 
struction  of  boundaries  (more  generally,  singularities)  in  the  object  function.  One  of  the 
sufficient  conditions  suggests  that  all  boundaries  in  the  object  function  can  stably  be  recon¬ 
structed  from  data  acquired  over  [0,  37r/2).  Also,  as  discussed  in  Appendix  A,  depending 
upon  the  boundary  locations  within  an  object  function,  it  is  theoretically  possible  that  all 
boundaries  can  stably  be  recovered  from  short-scan  or  7r-scheme  data.  More  important,  as 
the  results  of  numerical  studies  below  show,  under  realistic  practical  conditions  (e.g.,  in  the 
presence  of  data  noise),  images  reconstructed  from  short-scan  or  7T-scheme  data  generally 
appear  to  be  comparable  to  those  obtained  from  full-scan  data. 
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4  Results 


Using  simulation  studies,  we  numerically  investigated  the  accuracy  of  images  reconstructed 
from  reduced-scan  data  in  reflectivity  tomography. 


4.1  Data  and  Reconstruction  Algorithm 

We  used  three  different  numerical  phantoms,  which  are  shown  in  Fig.  3,  to  generate 
simulated  measurement  data  corresponding  to  the  imaging  transform  in  Eqn.  1.  The 
phantom  in  the  left  panel  of  Fig.  3  is  comprised  of  three  2D  Gaussian  functions,  representing 
a  smooth  object  function  containing  no  boundaries,  whereas  the  phantoms  in  the  middle 
and  right  panels  of  Fig.  3  contain  boundaries  and  complex  structures.  We  refer  to  the 
two  phantoms  in  the  left  and  middle  panels  of  Fig.  3  as  the  Gaussian  phantom  and  the 
ellipse  phantom,  respectively.  The  phantom  on  the  right  panel,  known  as  the  Shepp- 
Logan  phantom,  is  widely  used  in  the  medical  imaging  community.  When  we  simulated 
the  measured  data,  the  phantoms  were  centered  at  the  origin  of  the  x-y  coordinate  system, 
which  was  also  the  center  of  the  circular  source-receiver  trajectory  (see  Fig.  1).  The  matrix 
size  of  the  reconstructed  images  is  128x128. 

It  remains  unclear  whether  analytic  algorithms  exist  for  accurate  reconstruction  of  im¬ 
ages  from  reduced-scan  data.  In  this  work,  we  use  the  EM  algorithm  [11]  for  image  recon¬ 
struction.  For  notational  convenience,  let  x  =  (r,  6)  and  y  =  (£,  <j>)  denote  2D  vectors  in 
the  image  and  data  spaces,  respectively.  Equation  1  can  then  be  expressed  as 

g(y)=  [  dx  h(x,  y)  f(x)  MyeDy,  (14) 

J  D  x 


where  the  real  and  non-negative  functions  g{y)  and  f(x)  denote  the  data  and  object  func¬ 
tions  that  are  supported  on  the  domains  Dy  and  Dx,  respectively.  (For  instance,  for  a 
short-scan  configuration,  Dy  =  {y\  —  R  <  £  <  R,  $  6  7r}.)  The  kernel  h(x,y)  of  the 
imaging  transformation  is  given  by 


h(x,  y)  —  5  ^  r2  +  R2  —  2rRcos((f>  —  0)] 


We  use 

.  f^lx)  r  h(x,y)g(y )  .  . 

Z)  =  JD,dvh(x,y)  Jd,  V  SDJxh(x,y)fW(x)  (16) 

as  the  estimate  of  the  object  function,  where  n  is  the  number  of  iterations.  The  initial 
image  f^{x,y)  is  often  chosen  as  a  non-zero  constant  function.  The  EM  algorithm  has 
been  shown  to  yield  the  maximization-likelihood  solution  to  Eqn.  14  when  the  data  function 
g(y)  contain  Poisson  noise  [15].  From  a  practical  point  of  view,  the  EM  algorithm  is  easy 
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to  implement  because  it  involves  only  forward  and  backward  transformations  between  the 
image  and  data  spaces.  In  our  work,  in  an  attempt  to  enhance  the  computational  efficiency, 
we  utilized  an  ordered-subsets  [16]  implementation  of  the  EM  algorithm. 

4.2  Reconstruction  of  Smooth  Object  Functions 

We  first  investigated  the  accuracy  of  reconstructed  images  of  a  smooth  object  function  con¬ 
taining  no  boundaries.  To  do  so,  we  used  the  Gaussian  phantom  in  Fig.  3  to  generate  three 
noiseless  data  sets  containing  measurements  over  the  angular  ranges  of  [0,  7t/2),  [0,  37t/4), 
and  [0, 7r).  The  radius  of  the  source-receiver  trajectory  was  taken  to  be  R  =  r0,  where  r0 
represents  half  of  the  image-array  size.  From  the  three  data  sets,  we  reconstructed  images 
and,  in  each  case,  calculated  a  difference  image  between  the  reconstructed  image  and  the 
Gaussian  phantom. 

In  Fig.  4,  we  display  profiles  through  the  three  difference  images  at  x  =  64  and  y  =  44 
obtained  from  the  [0,  7t/2),  [0,  37t/4),  and  [0, 7r)  data  sets,  respectively.  It  can  be  observed 
that  the  difference  profiles  obtained  from  the  [0,  n)  data  set  are  almost  zero.  Because 
the  Gaussian  phantom  studied  has  no  obvious  symmetries,  such  an  observation  appears 
to  suggest  that  an  accurate  image  of  a  smooth  object  function  can  be  reconstructed  from 
short-scan  data.4  It  can  also  be  seen  in  Fig.  4  that  the  difference  profiles  calculated  from 
the  [0,  7t/2)  and  [0,  37t/4)  data  sets  are  significantly  different  from  zero,  suggesting  that 
accurate  images  cannot  be  reconstructed  from  data  acquired  over  angular  intervals  less 
than  7 r. 


4.3  Reconstruction  of  Images  Containing  Boundaries 

For  the  ellipse  and  Shepp-Logan  phantoms  that  contain  boundaries,  we  generated  simu¬ 
lated  data  sets  containing  measurements  over  angular  ranges  [0, 7r/2),  [k,  27t),  [7t/3,  27t/3]  U 
[7r, 47t/3]  U  [57t/3,  27 r],  [0,37t/2),  and  [0, 27r),  as  shown  in  Fig.  5.  The  second  and  third 
configurations  correspond  to  short-scan  and  7r-scheme  configurations,  whereas  the  last  con¬ 
figuration  is  the  full-scan  configuration.  For  each  of  the  ellipse  and  Shepp-Logan  phantoms 
and  for  each  of  the  scanning  configurations  in  Fig.  5,  we  generated  two  sets  of  noiseless  data 
by  use  of  two  circular  source-receiver  trajectories  with  different  radii.  Using  these  noiseless 
data  as  the  means,  we  generated  the  corresponding  noisy  data  by  adding  Gaussian  noise. 

In  Fig.  6,  we  display  images  reconstructed  from  noiseless  and  noisy  data  acquired  over 
[0, 7r/2)  (see  the  configuration  in  the  far  left  panel  of  Fig.  5).  Clearly,  all  of  the  images 
contain  significant  distortions,  suggesting  that  images  cannot  be  reconstructed  accurately 

4In  reflectivity  tomography,  the  stable  reconstruction  of  smooth  components  of  an  object  function  can 
be  deduced  by  exploiting  a  correspondence  between  the  Sobolev  wavefront  sets  of  the  data  and  object 
functions  [17,18]. 
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from  data  acquired  over  an  angular  interval  less  than  7r.  This  observation  is  consistent  with 
that  for  smooth  images  above. 

For  the  ellipse  phantom,  we  show  in  Figs.  7  and  8  images  reconstructed  from  data 
acquired  over  angular  ranges  [7r,  2n),  [7t/3, 27t/3]  U  [7r, 47t/3]  U  [57r/3,27r],  [0, 37r/2),  and 
[0,27t).  The  radii  of  the  source-receiver  trajectory  for  obtaining  the  results  in  Figs.  7  and 
8  were  r0  and  3r0,  respectively,  where  ro  represents  half  of  the  image-array  size.  Images 
reconstructed  from  noiseless  data  are  shown  in  the  first  columns  in  Figs.  7  and  8.  The 
overall  visual  quality  of  the  noiseless  images  reconstructed  from  short-scan  or  7r-scheme 
data  appears  comparable  to  that  obtained  from  full-scan  data.  In  particular,  for  the  case 
of  R  =  3r0  (Fig.  8),  which  may  correspond  to  a  typical  realistic  experimental  geometry, 
the  noiseless  images  reconstructed  from  the  short-scan  and  7r-scheme  data  appear  to  be 
virtually  identical  to  that  reconstructed  from  the  full-scan  data.  Images  in  the  second,  third, 
and  fourth  columns  in  Figs.  7  and  8  were'  reconstructed  from  data  containing  Gaussian 
noise  with  standard  deviations  of  1,  2,  and  4,  respectively.  These  noisy  images  appear 
distinct,  but  are  not  qualitatively  different  (for  a  given  noise  level).  Differences  in  the 
statistical  characteristics  of  the  images  (i.e.  the  appearance  of  the  image  noise)  are  to  be 
expected  because  the  same  reconstruction  algorithm  is  applied  to  different  noisy  data  sets. 
Furthermore,  the  full-scan  data  contain  redundant  information  that  is  not  present  in  the 
short-scan  or  7r-scheme  data,  which  serves  to  effectively  average  out  certain  components  of 
the  data  noise. 

For  the  Shepp-Logan  phantom,  we  show  images  in  Figs.  9  and  10  reconstructed  from 
data  acquired  over  angular  ranges  [n,  27t),  [7t/3,  27t/3]  U  [7t,  4tt/3]  U  [57r/3,  27t],  [0, 37r/2),  and 
[0,  27t),  respectively.  The  radii  of  the  source-receiver  trajectory  for  obtaining  the  results  in 
Figs.  9  and  10  were  r0  and  3 r0,  respectively.  Images  reconstructed  from  noiseless  data  are 
shown  in  the  first  columns  of  Figs.  9  and  10,  and  images  in  the  second,  third,  and  fourth 
columns  in  Figs.  9  and  10  were  reconstructed  from  data  containing  Gaussian  noise  with 
standard  deviations  of  1,  2,  and  4,  respectively.  Again,  for  images  in  Figs.  9  and  10,  one 
can  make  observations  similar  to  those  for  images  in  Figs.  7  and  8.  The  overall  visual 
quality  of  the  images  reconstructed  from  short-scan  or  7r-scheme  data  appears  comparable 
to  that  obtained  from  the  full-scan  data,  and  images  reconstructed  from  various  noisy  data 
sets  appear  qualitatively  similar  (for  a  given  noise  level). 

According  to  microlocal  analysis  [10],  certain  image  boundaries  may  not  theoretically  be 
reconstructed  stably  from  short-scan  or  7r-scheme  data.  For  example,  in  the  short-scan  case 
in  which  the  source-receiver  trajectory  is  within  the  lower  half  plane  (see  the  configuration 
in  the  second  panel  in  Fig.  5),  the  microlocal  analysis  suggests  that  vertical  boundaries 
in  the  ellipse  phantom  that  reside  in  the  upper  half  plane  cannot  be  reconstructed  stably. 
This  is  because  the  normal  vectors  of  such  vertical  boundaries  do  not  intersect  the  source- 
receiver  trajectory  in  the  lower  half  plane.  From  the  noiseless  image  in  the  far  left  panel  of 
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the  first  row  in  Fig.  7,  one  can  see  artifacts  near  the  vertical  boundaries  of  the  square  region 
in  the  upper  half.  However,  it  is  interesting  to  observe  that  such  artifacts  do  not  appear  to 
be  prominent.  Furthermore,  as  the  radius  of  the  source-receiver  trajectory  increases,  such 
artifacts  become  almost  undetectable  numerically,  as  confirmed  by  the  virtually  artifact- 
free  noiseless  images  in  the  far  left  panels  of  the  first  and  second  rows  in  Fig.  8  for  the  case 
R  —  3r0. 

The  boundaries  in  the  images  in  the  far  left  panels  of  the  third  rows  in  Figs.  7-10  appear 
identical  to  those  displayed  in  the  far  left  panels  of  the  fourth  rows  in  Figs.  7-10.  This 
confirms  our  assertion  (see  Appendix  A)  that  image  boundaries  can  stably  be  reconstructed 
from  data  acquired  over  an  angular  range  of  37t/2.  As  images  in  the  second,  third,  and 
fourth  columns  in  Figs.  7-10  show,  for  a  given  noise  level,  image  boundaries  reconstructed 
from  noisy  data  acquired  with  the  reduced-scan  and  full-scan  data  appear  qualitatively 
similar. 

For  completeness,  we  also  used  Norton’s  algorithm  (Eqns.  2,  3,  and  4)  to  reconstruct 
images  of  the  Shepp-Logan  phantom  from  data  that  are  acquired  with  the  short-scan  and 
full-scan  configurations  in  which  the  radius  of  the  source-receiver  trajectory  is  3ro.  Images 
reconstructed  from  short-scan  and  full-scan  data  are  shown  in  the  first  and  second  rows  in 
Fig.  11,  respectively.  Noiseless  images  are  shown  in  the  first  columns  in  Fig.  11,  whereas 
noisy  images  shown  in  the  second,  third,  and  fourth  columns  in  Fig.  11  were  reconstructed 
from  data  containing  Gaussian  noise  with  standard  deviations  of  1,  2,  and  5,  respectively. 
As  expected,  images  reconstructed  from  short-scan  data  contain  significant  artifacts.  Also, 
images  reconstructed  from  full-scan  data  contain  slight  ringing  artifacts  and  have  poorer 
resolution  as  compared  to  those  obtained  with  the  EM  algorithm,  suggesting  that  Norton’s 
algorithm  is  susceptible  to  numerical  errors. 

5  Discussion 

In  this  work,  we  used  a  potato-peeler  perspective  to  investigate  data  symmetries  in  reflectiv¬ 
ity  tomography.  Based  upon  the  identified  data  symmetries,  we  showed  heuristically  that 
an  object  function  with  a  compact  support  can  be  specified  by  knowledge  of  half  of  the  full- 
scan  data.  This  observation  has  led  to  the  development  of  reduced-scan  configurations  such 
as  short-scan  and  7r-scheme  configurations  in  reflectivity  tomography.  Reduced-scan  reflec¬ 
tivity  tomography  not  only  poses  a  theoretically  interesting  reconstruction  problem,  but 
also  has  practical  implications.  For  example,  in  certain  situations,  it  may  not  be  experimen¬ 
tally  possible  to  collect  data  over  a  full  angular  range  of  27r,  thus  demanding  reduced-scans. 
A  unique  feature  of  the  potato-peeler  perspective  is  that  it  provides  informative  insights 
into  a  complex  mathematical  problem  in  a  conceptually  straightforward  manner.  We  also 
employed  the  results  of  microlocal  analysis  to  analyze  theoretically  the  stability  of  recon- 
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structing  image  boundaries  from  reduced-scan  data,  and  we  derived  theoretically  sufficient 
conditions  on  reduced-scan  configurations  which  acquire  data  for  accurate  reconstruction 
of  an  object  function. 

As  predicted  by  microlocal  analysis,  certain  image  boundaries  cannot  stably  be  recon¬ 
structed  from  short-scan  or  7r-scheme  data  in  reflectivity  tomography.  Although  artifacts 
due  to  the  unstable  boundary  reconstruction  can  be  observed,  such  artifacts  appear  gen¬ 
erally  weak  and  become  difficult  to  discern  when  the  data  contain  noise  and/or  when  the 
radius  R  of  the  source-receiver  trajectory  becomes  large.  For  example,  for  R  =  3ro,  such 
artifacts  were  virtually  undetectable  for  the  ellipse  and  Shepp-Logan  phantoms  studied. 
This  is  significant  because  the  ellipse  phantom  containing  the  square  structures  was  specif¬ 
ically  designed  to  reveal  the  possible  instabilities  of  boundary  reconstruction  in  short-scan 
reflectivity  tomography  predicted  by  the  microlocal  analysis,  and  because  the  Shepp-Logan 
phantom  is  sufficiently  complex  to  represent  many  real-world  reflecting  objects  that  may 
be  of  practical  interest.  Conceptually,  these  results  might  be  explained  by  the  fact  that 
the  data  are  becoming  more  “Radon-like”  as  the  radius  R  of  the  source-receiver  trajectory 
increases.  (In  the  limit  R  — >•  oo,  Eqn.  1  reduces  to  the  Radon  transform,  and  an  exact 
image  can  be  reconstructed  from  short-scan  or  7r-scheme  data.)  Our  numerical  results  have 
important  practical  implications:  For  many  reflecting  objects  possessing  boundaries  and 
for  scanning  configurations  with  different  R  that  may  arise  in  realistic  experiments,  images 
reconstructed  from  reduced-scan  data,  such  as  the  short-scan  and  7r-scheme  data,  can  have 
a  numerical  accuracy  similar  to  that  of  full-scan  images.  The  numerical  results  also  ver¬ 
ified  the  sufficient  condition  that  all  image  boundaries  can  stably  be  reconstructed  from 
reduced-scan  data  taken  over  an  angular  range  of  37t/2.  For  smooth  object  functions,  our 
numerical  results  suggest  that  images  can  be  reconstructed  accurately  from  short-scan  and 
7r-scheme  data. 

We  are  currently  developing  a  mathematical  formulation  for  the  potato-peeler  perspec¬ 
tive  for  reflectivity  tomography  and  will  report  such  results  elsewhere.  The  perspectives 
and  methods  utilized  in  this  work  can  readily  be  extended  to  investigation  of,  e.g.,  recon¬ 
struction  problems  in  reflectivity  tomography  for  which  the  source-receiver  trajectory  is 
not  a  circular  curve,  but  is  a  general  curve  satisfying  certain  conditions.  This  work  is  also 
relevant  to  generalized  Radon  transforms  that  integrate  over  certain  (non-circular)  curves 
and  to  exterior  reconstruction  problems.  Investigation  of  these  topics  is  currently  under 
way. 
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Appendix  A:  Stability  of  Reconstructing  Image  Bound¬ 
aries  from  Reduced-Scan  Data 

Microlocal  analysis  [14,18,19]  can  be  employed  for  examining  whether  image  boundaries 
can  stably  be  reconstructed  from  reduced-scan  data  in  reflectivity  tomography.  Without 
loss  of  generality,  we  consider  boundary  reconstruction  from  short-scan  data  acquired  with 
the  configurations  shown  in  Fig.  12,  where  the  source-receiver  trajectories  are  semi-circles 
of  radii  P  residing  in  the  lower  half  plane.  According  to  the  microlocal  analysis  [10],  a 
boundary  can  stably  be  reconstructed  if  and  only  if  the  normal  vector  at  that  boundary 
intersects  the  source-receiver  trajectory  [17].  For  example,  as  shown  in  Fig.  12a,  assume 
R  =  Ri,  and  consider  the  three  boundaries  labeled  Pi,  P2 ,  and  P3.  The  normal  vector 
at  boundary  Pi  intersects  the  source-receiver  trajectory,  and  therefore  boundary  Pi  can 
stably  be  reconstructed.  However,  because  the  normal  vectors  at  boundaries  P2  and  P3 
do  not  intersect  the  trajectory  with  a  radius  Pi,  boundaries  P2  and  P3  cannot  stably 
be  reconstructed.  If  the  radius  of  the  source-receiver  trajectory  is  increased  to  be  large 
enough,  the  normal  vectors  that  are  not  parallel  to  the  x-axis  will  eventually  intersect  the 
trajectory,  and  the  number  of  boundaries  that  can  stably  be  reconstructed  will  increase.  For 
example,  as  shown  in  Fig.  12,  when  P  =  P2,  boundary  P2  can  now  be  reconstructed  stably. 
Therefore,  when  one  reconstructs  an  image  from  short-scan  data  in  reflectivity  tomography, 
the  manifestation  of  artifacts  due  to  boundaries  that  cannot  be  stably  reconstructed  would 
depend  on  the  magnitude  of  P.  Notice  that  the  normal  vector  at  boundary  P3  is  parallel 
to  the  rr-axis,  and  will  therefore  never  intersect  the  source-receiver  trajectory  with  a  finite 
radius  in  the  lower  half  plane.  Therefore,  one  can  conclude  that  certain  boundaries  cannot, 
in  theory,  be  reconstructed  stably  from  short-scan  data  in  reflectivity  tomography.  Similar 
conclusions  can  be  made  for  the  stability  of  boundary  reconstruction  in  7r-scheme  reflectivity 
tomography.  However,  as  the  results  in  numerical  studies  indicated,  for  object  functions 
with  complicated  boundaries  and  for  source-receiver  trajectories  with  radii  P  of  practical 
interest,  artifacts  due  to  unstable  reconstruction  of  the  boundaries  from  short-scan  and 
7r-scheme  data  do  not  appear  to  be  prominent. 

Based  upon  the  microlocal  analysis,  one  can  derive  theoretically  sufficient  conditions 
for  stable  reconstruction  of  image  boundaries.  In  the  limiting  case,  i.e.,  P  ->  00,  the  data 
function  (see  Eqn.  1)  in  reflectivity  tomography  reduces  to  the  Radon  transform,  and 
normal  vectors  of  any  boundaries  will  intersect  the  source-receiver  trajectory  at  infinity. 
Therefore,  all  of  the  boundaries  can  stably  be  reconstructed  from  short-scan  or  7r-scheme 
Radon  data.  Also,  considering  the  scanning  configuration  in  the  fourth  panel  of  Fig.  5, 
the  source-receiver  trajectory  of  radius  P  in  this  configuration  covers  an  angular  interval  of 
37t/2.  In  this  case,  it  can  readily  be  verified  that  any  normal  vectors  on  image  boundaries 
intersect  the  source-receiver  trajectory.  Therefore,  one  can  obtain  a  sufficient  condition 
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that  all  of  the  boundaries  can  stably  be  reconstructed  from  data  acquired  over  a  reduced 
angular  interval  that  spans  37r/2.  Furthermore,  one  can  design  a  scanning  configuration 
such  that  all  of  the  boundaries  within  an  object  function  can  stably  be  reconstructed  from, 
e.g.,  the  short-scan  data.  As  shown  in  Fig.  12b,  for  the  entire  object  that  is  placed  in 
the  lower  half  of  the  plane  and  is  enclosed  by  the  source-receiver  trajectory,  all  boundaries 
in  the  object  can  stably  be  reconstructed  because  their  normal  vectors,  in  this  situation, 
intersect  the  source-receiver  trajectory. 
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Figure  Captions 

1.  Schematic  illustration  of  the  scanning  geometry  in  reflectivity  tomography. 


2.  Schematic  illustration  of  the  measurements  at  the  edge  of  an  object  function  in  reflec¬ 
tivity  tomography  (left  panel).  Points  A,  B,  and  C  in  the  object  function  contribute  to  the 
measured  data  on  curves  with  the  same  labels  in  the  data  space  (right  panel). 

3.  Numerical  phantoms  used  in  simulation  studies.  The  phantom  in  the  left  panel,  which  is 
referred  to  as  the  Gaussian  phantom,  is  comprised  of  three  Gaussian  functions,  representing 
a  smooth  object  function.  The  phantoms  in  the  middle  and  right  panels,  which  contain 
boundaries  (i.e.,  discontinuities),  are  referred  to  as  the  ellipse  and  Shepp-Logan  phantoms, 
respectively,  and  are  used  to  investigate  boundary-reconstruction  stability  from  short-scan 
and  7r-scheme  data  in  reflectivity  tomography. 

4.  Profiles  through  difference  images  at  x  =  64  (left)  and  y  =  44  (right)  obtained  from 
data  acquired  over  [0,  tt)  (solid),  [0, 37r/4)  (dotted),  and  [0,7r/2)  (dashed). 

5.  Scanning  configurations  that  acquire  data  over  angular  intervals  [0, 7r/2)  (left),  [tt,  27t) 
(second  from  left),  [7r/3, 27r/3]U[7r,  47r/3]u[57r/3,  2tv]  (middle),  [0,  37t/2)  (second  from  right), 
and  [0, 27 r)  (right).  The  configurations  in  the  second,  third,  and  fifth  panels  are  also  referred 
to  as  the  short-scan,  7r-scheme,  and  full-scan  configurations,  respectively. 

6.  Images  reconstructed  from  noiseless  and  noisy  data  acquired  over  [0,  7t/2)  for  the  ellipse 
and  Shepp-Logan  phantoms.  The  radii  of  the  source-receiver  trajectory  were  set  to  r0  in  (a) 
and  (c),  and  3 r0  in  (b)  and  (d).  Images  in  the  first  column  were  reconstructed  from  noiseless 
data,  and  images  in  the  second,  third,  and  fourth  columns  were  reconstructed  from  data 
containing  three  different  levels  of  Gaussian  noise.  All  of  these  images  contain  significant 
distortions,  suggesting  that  images  cannot  accurately  be  reconstructed  from  data  acquired 
over  angular  intervals  less  than  tt. 

7.  Images  reconstructed  for  the  ellipse  phantom  from  data  acquired  over  angular  intervals 
[%,  2tt)  (1st  row),  [7r/3, 27r/3]  U  [7t,47t/3]  U  [57t/3,  27t]  (2nd  row),  [0,  37t/2)  (3rd  row),  and 
[0,  27t)  (4th  row).  The  radius  of  the  source-receiver  trajectory  is  r0.  Images  in  the  first 
column  were  reconstructed  from  noiseless  data,  and  images  in  the  second,  third,  and  fourth 
columns  were  reconstructed  from  data  containing  three  different  levels  of  Gaussian  noise. 

8.  Images  reconstructed  for  the  ellipse  phantom  from  data  acquired  over  angular  intervals 
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[tt,  27t)  (1st  row),  [7r/3, 27r/3]  U  [7r,  47r/3]  U  [57r/3,27r]  (2nd  row),  [0,  37t/2)  (3rd  row),  and 
[0,  2tt)  (4th  row).  The  radius  of  the  source-receiver  trajectory  is  3r0.  Images  in  the  first 
column  were  reconstructed  from  noiseless  data,  and  images  in  the  second,  third,  and  fourth 
columns  were  reconstructed  from  data  containing  three  different  levels  of  Gaussian  noise. 

9.  Images  reconstructed  for  the  Shepp-Logan  phantom  from  data  acquired  over  angular 
intervals  [7t,27t)  (1st  row),  [7r/3,  27t/3] U [7t,  47t/3]  U [57t/3,  27t]  (2nd  row),  [0,  37t/2)  (3rd  row), 
and  [0, 2n)  (4th  row).  The  radius  of  the  source-receiver  trajectory  is  r0.  Images  in  the  first 
column  were  reconstructed  from  noiseless  data,  and  images  in  the  second,  third,  and  fourth 
columns  were  reconstructed  from  data  containing  three  different  levels  of  Gaussian  noise. 

10.  Images  reconstructed  for  the  Shepp-Logan  phantom  from  data  acquired  over  angular 
intervals  [7r,  27r)  (1st  row),  [7r/3,  27t/3]  U [7t,  47t/3]  U [57t/3,  27t]  (2nd  row),  [0,  37t/2)  (3rd  row), 
and  [0,27t)  (4th  row).  The  radius  of  the  source-receiver  trajectory  is  3r0.  Images  in  the 
first  column  were  reconstructed  from  noiseless  data,  and  images  in  the  second,  third,  and 
fourth  columns  were  reconstructed  from  data  containing  three  different  levels  of  Gaussian 
noise. 

11.  Images  reconstructed  by  use  of  Norton’s  algorithm  from  short-scan  data  (1st  row) 
and  full-scan  data  (2nd  row)  that  were  acquired  with  a  circular  source-trajectory  of  radius 
R  =  3r0.  Images  in  the  first  column  were  reconstructed  from  noiseless  data,  and  images  in 
the  second,  third,  fourth  columns  were  reconstructed  from  data  containing  three  different 
levels  of  Gaussian  noise. 

12.  Short-scan  configurations  and  boundaries  in  the  objects,  (a)  Some  of  the  boundaries 
such  as  P2  and  P3  reside  in  the  upper  half  plane;  and  (b)  the  entire  object  resides  in  the 
lower  half  plane  and  is  enclosed  by  the  source-receiver  trajectory. 
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